有限差分法和托马斯算法（或追赶法）对一个二阶线性边值问题进行数值求解

2744557306

使用有限差分法和托马斯算法（或追赶法）对一个二阶线性边值问题进行数值求解。这种方法通常用于数值解微分方程。
% \* m+ q2 Q( }6 q" {& W以下是代码的简要解释：

1.使用 inline 函数定义了三个函数 p(x)、q(x) 和 r(x)，它们表示微分方程的系数。
9 \" ]1 @4 y" R* K; S7 j2.设置了参数，如间隔数 N、初始和边界条件 a0、b0、af、bt 以及间隔大小 h。
3.基于微分方程的有限差分离散化，计算了系数 a、b、c 和 d。
4.使用托马斯算法（或追赶法）解决了三对角方程组。
5.将结果与由数组 zj 表示的解析解进行了比较。
6.将数值解和解析解并排显示，以便比较。

1. p=inline('-2/x');8 l2 a/ Y& y2 ?9 P9 u
   
2.   q=inline('2/x^2');7 m. O0 l7 z5 i9 T, L, Q7 }! k
   
3.   r=inline('sin(log10(x)/log10(exp(1)))/x^2');
   ) t4 V5 k, D+ z, V
4.   N=9;
   5 H4 P: U# z\" A! p
5.   a0=1;b0=2;& V/ G0 d8 S+ F
   
6.   af=1;bt=2;
     n6 h- P7 M\" y. T& k. G
7.   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%9 ~1 v2 q4 s, ^& d! S. u6 D
   
8.   h=(b0-a0)/(N+1);
   * N& M( E6 G  {3 H
9.   x=a0+h;
   8 ^1 Y% X) T) P2 D+ x( N7 P, E
10.   a(1)=2+h*h*q(x);7 n* `\" {3 z, q. L
    
11.   b(1)=-1+(h/2)*p(x);
    , A7 \8 C8 l/ U7 K% ]8 q' P5 z
12.   d(1)=-h*h*r(x)+(1+(h/2)*p(x))*af;- \* ?' ?$ r$ v* \* C
    
13.     for i=2:N-1
    ( ]  @5 ?7 ~7 j
14.         x=a0+i*h;7 O. p/ X! B1 I) {+ n1 U5 P* C# Q
    
15.         a(i)=2+h*h*q(x);
    $ B9 f) ~; F/ U
16.         b(i)=-1+(h/2)*p(x);/ e( _$ X' I# M/ q' b
    
17.         c(i)=-1-(h/2)*p(x);
    / l: c* j9 G2 x. ~# F
18.         d(i)=-h*h*r(x);
    & ]) S0 b6 B! ?9 |7 O' W# g
19.      end
    7 k+ D, M  S! z. C' u\" p
20.      x=b0-h;4 g! b, n. M% h7 I( v) @
    
21.      a(N)=2+h*h*q(x);
    1 ]# q, Z- g! k( o6 l
22.      c(N)=-1-(h/2)*p(x);
    / {5 l8 l: _/ n
23.      d(N)=-h*h*r(x)+(1-(h/2)*p(x))*bt;  `* M7 w7 a2 Z  l
    
24.      %%%%%%%%%追赶法%%%%%%%%%%%%%%%%%%
    4 Z, b, v# [9 }: d+ u0 S9 H
25.      %y=trisys(c,a,b,d)
    ' T. M  C# D  ~! s9 A1 M1 K
26.      L(1)=a(1);; v7 K\" L- q, v! z4 c. u2 X* B( w6 ^
    
27.      u(1)=b(1)/a(1);
    \" X/ x2 c$ c, i/ J
28.      for i=2:N-1
    ; _2 W& Q! Z! c: R5 h8 M
29.          L(i)=a(i)-c(i)*u(i-1);8 o* V9 O  t  N$ _
    
30.          u(i)=b(i)/L(i);0 X9 m% h2 r) g. g. s1 z, V
    
31.      end) X+ C( P( x1 P5 P; o* v
    
32.      L(N)=a(N)-c(N)*u(N-1);) l3 N' K' E: Y0 ^2 {) X; s
    
33.      z(1)=d(1)/L(1);
    % g* @% ]* E; h
34.      for i=2:N
    ! p& S/ X* W2 K
35.          z(i)=(d(i)-c(i)*z(i-1))/L(i);
    7 r' b( ~; x  X3 e
36.      end
    : ?' `3 d/ b; b1 p* S7 q3 V5 N6 C
37.      y(N)=z(N);' Y9 d, u3 r* l! }0 u
    
38.      for i=N-1:-1:1: c1 U3 J# @1 w* D, P  k
    
39.          y(i)=z(i)-u(i)*y(i+1);$ r/ L- _9 e. l0 M+ ^
    
40.      end8 |$ @2 c2 H1 k# a
    
41.      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    6 @; F7 J& n5 I\" f
42.      Y=[af,y,bt];
    : P( h1 o3 W+ V! G5 \
43. for i=1:N+2
    $ I, r! m  x$ j
44.       x=a0+(i-1)*h;
    ; ]6 C9 r7 a5 z0 `# o
45.       zj(i)=1.1392070132*x-0.03920701320/x^2-3*sin(log10(x)/log10(exp(1)))/10-cos(log10(x)/log10(exp(1)))/10;: D5 `! I7 g! X
    
46. end: G, B& A$ r) D2 n) k; U
    
47. disp('下面两列分别是数值解和近似解');
    + F3 h& K) K! Y: N+ X
48. re=[Y'      zj']

复制代码

/ Z: s& \6 ^5 i) x6 ^