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

2744557306

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

" \% o8 C5 W& W, o- t7 f1.使用 inline 函数定义了三个函数 p(x)、q(x) 和 r(x)，它们表示微分方程的系数。
2.设置了参数，如间隔数 N、初始和边界条件 a0、b0、af、bt 以及间隔大小 h。
6 a( I$ t- K, O. l3 r( y3.基于微分方程的有限差分离散化，计算了系数 a、b、c 和 d。
' d% O# X, b7 S9 J1 m6 X9 _, B1 i% k4.使用托马斯算法（或追赶法）解决了三对角方程组。
5.将结果与由数组 zj 表示的解析解进行了比较。
# {2 u9 ?9 D8 f2 k: Z6.将数值解和解析解并排显示，以便比较。

1. p=inline('-2/x');/ L. L$ c3 i( K\" y
   
2.   q=inline('2/x^2');* [! Q! P$ C8 }. J) J4 J8 s/ Y
   
3.   r=inline('sin(log10(x)/log10(exp(1)))/x^2');
   + q% V; a5 g; Q% H
4.   N=9;
   # h3 X8 o! t( j
5.   a0=1;b0=2;
   # R6 a2 |6 n* g7 M# O
6.   af=1;bt=2;) B* w8 s& p6 c( _2 M
   
7.   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   7 Z4 L/ {! F* c# @0 R
8.   h=(b0-a0)/(N+1);6 s' W- o( t3 `$ ?9 f- z
   
9.   x=a0+h;
   3 E7 Y$ b3 K+ Q, t  Q
10.   a(1)=2+h*h*q(x);\" ^! i1 ]/ v1 w/ ^
    
11.   b(1)=-1+(h/2)*p(x);* N& {1 M: x$ ?\" A7 p\" }( }
    
12.   d(1)=-h*h*r(x)+(1+(h/2)*p(x))*af;' H, F1 V. v; t. V; ~\" H1 j
    
13.     for i=2:N-1
    8 q1 ?' y8 J* Z: [\" F, b8 D. _. @
14.         x=a0+i*h;
    , }. `7 t& v  r+ C  ]1 D
15.         a(i)=2+h*h*q(x);
    6 s9 k$ f  ~+ Q  R* P+ h7 ?
16.         b(i)=-1+(h/2)*p(x);7 N) b: G8 k' ]0 i
    
17.         c(i)=-1-(h/2)*p(x);9 J) f1 N9 L3 C9 F: X) k+ I1 ~
    
18.         d(i)=-h*h*r(x);
    2 [! `7 K7 x  s# t; w+ S7 c4 I
19.      end
    ) h\" l9 z) Z' H( p5 k/ o9 f
20.      x=b0-h;* X( t: v. B- h4 F5 H  j* [. T2 b
    
21.      a(N)=2+h*h*q(x);
    * J. n2 M6 r& f9 B4 y' @; ^
22.      c(N)=-1-(h/2)*p(x);3 P1 i; z8 J2 Q
    
23.      d(N)=-h*h*r(x)+(1-(h/2)*p(x))*bt;
    ; B% b/ @# b7 l1 g2 `
24.      %%%%%%%%%追赶法%%%%%%%%%%%%%%%%%%3 ~: h9 c0 H/ U* e
    
25.      %y=trisys(c,a,b,d)
    8 W5 p+ H\" y+ ]7 W! Q; i
26.      L(1)=a(1);
    0 w# R( g* x\" n$ X' ?
27.      u(1)=b(1)/a(1);$ I* H7 s3 a; Z# T, C/ F
    
28.      for i=2:N-1
    6 F8 @9 g/ `) J, x! s\" |$ f6 @
29.          L(i)=a(i)-c(i)*u(i-1);
    ' m/ C3 V2 t8 h' E0 _1 P, s
30.          u(i)=b(i)/L(i);0 q4 n# k1 O) m/ s: Z
    
31.      end
    8 N8 d0 G! F& w+ |2 W3 Z9 Q7 r2 t
32.      L(N)=a(N)-c(N)*u(N-1);: i6 p; g7 q$ C; e- u) C2 v
    
33.      z(1)=d(1)/L(1);
    ; T9 Q- z, g( ?
34.      for i=2:N& T; v, u3 p, R5 Q\" R: R, |/ w
    
35.          z(i)=(d(i)-c(i)*z(i-1))/L(i);
    9 [, {& |4 Q* E6 b# Q& I
36.      end8 X: N& }9 z; B
    
37.      y(N)=z(N);5 y7 ^* a/ i  y5 x\" E
    
38.      for i=N-1:-1:1
    9 K. f  R% V( Q
39.          y(i)=z(i)-u(i)*y(i+1);3 T: N6 O2 Y( d4 v
    
40.      end( |, V! w; u0 j8 Q9 B
    
41.      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    $ B' P5 G: F# _
42.      Y=[af,y,bt];% a1 r, R$ Q7 C\" e$ S
    
43. for i=1:N+2
    7 U0 e4 v( y3 D8 B. A( _; D
44.       x=a0+(i-1)*h;' y; Y& e; P- S9 B
    
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;, V0 B( V6 E5 v, z% ?3 f! N
    
46. end
    : b# X) l! T' e
47. disp('下面两列分别是数值解和近似解');
    , C6 c3 S\" `, f( P$ W8 R
48. re=[Y'      zj']

复制代码