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

2744557306

使用有限差分法和托马斯算法（或追赶法）对一个二阶线性边值问题进行数值求解。这种方法通常用于数值解微分方程。
2 l( _/ u7 ^' P& m* q! b" {  [以下是代码的简要解释：
& W4 R( |7 e6 U! ]
1.使用 inline 函数定义了三个函数 p(x)、q(x) 和 r(x)，它们表示微分方程的系数。
* F! x" a! O7 C( K8 |2.设置了参数，如间隔数 N、初始和边界条件 a0、b0、af、bt 以及间隔大小 h。
3.基于微分方程的有限差分离散化，计算了系数 a、b、c 和 d。
4.使用托马斯算法（或追赶法）解决了三对角方程组。
5.将结果与由数组 zj 表示的解析解进行了比较。
6.将数值解和解析解并排显示，以便比较。

1. p=inline('-2/x');
   \" Q' \3 P; b9 o! `( B
2.   q=inline('2/x^2');1 f  c8 I6 V. ]\" S( {
   
3.   r=inline('sin(log10(x)/log10(exp(1)))/x^2');
   ' x$ i. G- q. B4 L& y2 S0 n9 D7 _
4.   N=9;
   2 \& |$ z+ e  V; L) p
5.   a0=1;b0=2;& v# x! Q: L) Y) I: ^0 ^1 C
   
6.   af=1;bt=2;2 w! a9 A% z$ v3 k
   
7.   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%8 p1 N' l; A- [8 [$ ^; o  D
   
8.   h=(b0-a0)/(N+1);
   ' W4 z2 t* E# ?
9.   x=a0+h;9 `# ?- `2 y1 _( _# O& ]
   
10.   a(1)=2+h*h*q(x);/ {% d) k: E% A9 H7 `0 v
    
11.   b(1)=-1+(h/2)*p(x);5 g: _; p2 L* F% F\" n1 k
    
12.   d(1)=-h*h*r(x)+(1+(h/2)*p(x))*af;
    ' R& m2 a1 E\" ~
13.     for i=2:N-1! n. Z' C5 C  M: z  N
    
14.         x=a0+i*h;
      S; n) j! W( ]0 e( F7 B
15.         a(i)=2+h*h*q(x);
    ' G1 I0 P- \% j1 _% G
16.         b(i)=-1+(h/2)*p(x);7 N% K* g8 M/ T& n5 `2 V
    
17.         c(i)=-1-(h/2)*p(x);
    - u+ _\" x, n: f: g  q
18.         d(i)=-h*h*r(x);
      y/ g* k/ c5 L; [- l' ]\" Y. ]4 ]: P
19.      end1 E, L1 v8 g6 [. S% v2 G* A
    
20.      x=b0-h;
    9 p3 }5 T' `: a$ P: [; |
21.      a(N)=2+h*h*q(x);$ i; s: \: _% `* }( i* ^, f7 r! b
    
22.      c(N)=-1-(h/2)*p(x);
    5 r) F6 ^' F5 J0 _/ }
23.      d(N)=-h*h*r(x)+(1-(h/2)*p(x))*bt;$ P7 V2 \2 p% M2 ]
    
24.      %%%%%%%%%追赶法%%%%%%%%%%%%%%%%%%: J% y/ t+ l% S+ F
    
25.      %y=trisys(c,a,b,d)
    6 `  e' X. _. {+ @0 e9 k$ X
26.      L(1)=a(1);7 G3 l3 G, M0 \' F4 B6 U
    
27.      u(1)=b(1)/a(1);
    ! x  i' y' i2 K- l* ?8 N8 S! K4 z* ^
28.      for i=2:N-1
    : P9 e' Y% b# V\" e# W5 b
29.          L(i)=a(i)-c(i)*u(i-1);; R7 I: Y  t! B0 R\" J
    
30.          u(i)=b(i)/L(i);
    % n0 w5 J( P% A% L$ ~
31.      end
    1 y& p  \4 X) E3 p% u
32.      L(N)=a(N)-c(N)*u(N-1);, {% {0 e2 n$ T% Z& q% Y
    
33.      z(1)=d(1)/L(1);
    7 H# R6 g' P( }. S5 s
34.      for i=2:N
    - M- d; {- Q% j: u\" P
35.          z(i)=(d(i)-c(i)*z(i-1))/L(i);
    % d- E/ {9 b9 q\" \  b8 e
36.      end
    ; \( I/ U8 D# S' \
37.      y(N)=z(N);7 z9 ~; C  W$ @& t) ?( Z% D
    
38.      for i=N-1:-1:1& t. s* a9 d- S
    
39.          y(i)=z(i)-u(i)*y(i+1);1 ]1 l) a* L. W
    
40.      end- M& ?9 m5 k# {% T, W# g
    
41.      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    - l  G8 o# \5 Z; @5 c9 u) }
42.      Y=[af,y,bt];
    * ^5 ^1 b/ Y1 Q! b6 T1 v
43. for i=1:N+2
    \" i0 [4 Z8 B& r) O
44.       x=a0+(i-1)*h;
    & t/ N/ k! D# B6 x1 G
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;
    9 Y  {: c# }4 e) q
46. end
    4 g) }1 @2 J! l2 g( B% ^
47. disp('下面两列分别是数值解和近似解');
    , I& S' @% D. H% q
48. re=[Y'      zj']

复制代码

) ~8 L9 t# g: c+ r9 s, O3 C