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

2744557306

使用有限差分法和托马斯算法（或追赶法）对一个二阶线性边值问题进行数值求解。这种方法通常用于数值解微分方程。
以下是代码的简要解释：

9 U! r8 g6 j2 K& N) }3 E8 a1.使用 inline 函数定义了三个函数 p(x)、q(x) 和 r(x)，它们表示微分方程的系数。
9 M" O: Y+ p% c' ^2.设置了参数，如间隔数 N、初始和边界条件 a0、b0、af、bt 以及间隔大小 h。
3.基于微分方程的有限差分离散化，计算了系数 a、b、c 和 d。
1 a1 ?1 t2 J& O& _' Z; i" O4.使用托马斯算法（或追赶法）解决了三对角方程组。
5.将结果与由数组 zj 表示的解析解进行了比较。
% Y) V3 f/ |# a6 v6.将数值解和解析解并排显示，以便比较。

1. p=inline('-2/x');
   ) w+ Y\" Q% K$ f  n+ }) A# T: I6 d: x
2.   q=inline('2/x^2');* F9 Z  j3 B. m& i/ \$ f
   
3.   r=inline('sin(log10(x)/log10(exp(1)))/x^2');
   % M) ]  `/ C7 j
4.   N=9;
   / j( [$ D3 D+ r  A
5.   a0=1;b0=2;) Y\" r+ y9 R4 @4 [
   
6.   af=1;bt=2;
   % P4 M, U7 R% w: f
7.   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%' |' \2 ~3 y6 D: H: w) E
   
8.   h=(b0-a0)/(N+1);
   ; Z# q) |6 y\" a1 e2 H. V
9.   x=a0+h;' l  u$ i, ~/ V$ A! ?
   
10.   a(1)=2+h*h*q(x);$ L7 ~* {: C. P- k, q% Y
    
11.   b(1)=-1+(h/2)*p(x);- h  }3 X& c: r% ^$ F, U
    
12.   d(1)=-h*h*r(x)+(1+(h/2)*p(x))*af;9 t/ Z: @) n2 ?8 w
    
13.     for i=2:N-1+ V' \8 k: T% t$ r* `6 D/ U
    
14.         x=a0+i*h;
    3 O8 ]; p7 Z3 Y  _# X
15.         a(i)=2+h*h*q(x);
    + E0 y) G! U9 _& n
16.         b(i)=-1+(h/2)*p(x);
    * B) \# p) f2 D) X1 p/ U. W+ g; ~
17.         c(i)=-1-(h/2)*p(x);
    ) r! P5 V5 t7 |9 ]$ Q
18.         d(i)=-h*h*r(x);
    5 k* s/ N7 o6 u. S: ^
19.      end7 P5 c% ^- F' l% q
    
20.      x=b0-h;9 u* |1 ]1 p  |- ]1 [3 Y; g- {. [\" O5 {
    
21.      a(N)=2+h*h*q(x);7 u, I0 x1 T- i1 c1 Z0 ?2 @
    
22.      c(N)=-1-(h/2)*p(x);
    6 n+ _6 y' T; s( H- V2 f
23.      d(N)=-h*h*r(x)+(1-(h/2)*p(x))*bt;
    - x* N( A+ V) a. H9 r4 S- ^) `
24.      %%%%%%%%%追赶法%%%%%%%%%%%%%%%%%%
    ; Q6 f* _# v, v( c\" u( n
25.      %y=trisys(c,a,b,d)
    ( O5 m+ E* H8 A8 d! H! h4 a
26.      L(1)=a(1);
    \" y& q1 t% O2 g3 f
27.      u(1)=b(1)/a(1);- `2 t# \0 d, g* p, v/ W7 v0 J
    
28.      for i=2:N-1; U0 x# O0 l9 E: @# R  `
    
29.          L(i)=a(i)-c(i)*u(i-1);# m5 ]# f' O4 L5 X7 X
    
30.          u(i)=b(i)/L(i);
    / Q- U  W3 e! ]* s# O1 v. m
31.      end1 ~* R+ y! S/ m6 ^3 q- T
    
32.      L(N)=a(N)-c(N)*u(N-1);
    0 _/ G6 u. F  t$ v
33.      z(1)=d(1)/L(1);: l7 [\" v3 E\" T3 n1 B5 d
    
34.      for i=2:N8 q5 |: T+ t, @4 S
    
35.          z(i)=(d(i)-c(i)*z(i-1))/L(i);
    $ d\" G$ b( ]6 g# R* X
36.      end
    7 O& f6 y4 b: v0 f
37.      y(N)=z(N);
    + B/ i' F: _' j' K' w
38.      for i=N-1:-1:1
    * d5 R) E0 [/ g/ J
39.          y(i)=z(i)-u(i)*y(i+1);
    : N8 w& V9 W' ~& I8 W\" I$ p( }. H
40.      end- R0 _2 U( J) |+ l8 x7 v) c* E
    
41.      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    $ ]- J\" R: d2 C9 x9 ]: ]
42.      Y=[af,y,bt];* X8 e5 b, K* a2 q, J
    
43. for i=1:N+2( k0 |- X3 e: V4 A# r! j8 X4 p1 p
    
44.       x=a0+(i-1)*h;% O: ], G\" B; q\" N( ~
    
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;
    , |/ x# G7 K- F0 D. ]0 L/ q
46. end
    \" o: Q' I2 E) ]6 l; M5 E
47. disp('下面两列分别是数值解和近似解');
    : k$ w, T1 [0 z
48. re=[Y'      zj']

复制代码