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

2744557306

使用有限差分法和托马斯算法（或追赶法）对一个二阶线性边值问题进行数值求解。这种方法通常用于数值解微分方程。
以下是代码的简要解释：
7 [2 G$ c0 [5 r- x/ w0 f
1.使用 inline 函数定义了三个函数 p(x)、q(x) 和 r(x)，它们表示微分方程的系数。
, e2 T" |" z; k, s% R2.设置了参数，如间隔数 N、初始和边界条件 a0、b0、af、bt 以及间隔大小 h。
3.基于微分方程的有限差分离散化，计算了系数 a、b、c 和 d。
4.使用托马斯算法（或追赶法）解决了三对角方程组。
5.将结果与由数组 zj 表示的解析解进行了比较。
6.将数值解和解析解并排显示，以便比较。

1. p=inline('-2/x');\" Y- s: O: l- W, v; G  E
   
2.   q=inline('2/x^2');
   ; k! L\" ~' a3 R
3.   r=inline('sin(log10(x)/log10(exp(1)))/x^2');& t! g8 Z. ^% u2 g- v6 u& `( [
   
4.   N=9;
   * a  X; |3 K% F# V% Z1 j
5.   a0=1;b0=2;, Z- U/ B8 {# [8 w) d9 W2 e9 P
   
6.   af=1;bt=2;
     g7 u8 \; f: I; o
7.   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   * x6 z/ z6 _% |+ L* u: X) f5 Y; @; P
8.   h=(b0-a0)/(N+1);0 F. h5 s* Y- d/ x% N
   
9.   x=a0+h;* k, J# S2 ?, _+ o  R& D& I
   
10.   a(1)=2+h*h*q(x);
    , G4 s4 |' H* B5 L
11.   b(1)=-1+(h/2)*p(x);
    6 b2 D) W1 [0 K) R! ]$ V
12.   d(1)=-h*h*r(x)+(1+(h/2)*p(x))*af;
    9 o9 n0 {! n4 `/ K\" b: \
13.     for i=2:N-1
    8 T) `: ]$ J8 A+ N9 U% ^
14.         x=a0+i*h;6 ?/ w1 C% V% n# o
    
15.         a(i)=2+h*h*q(x);
    . }8 X; q( T( L$ g
16.         b(i)=-1+(h/2)*p(x);. \2 |5 A\" x, k  v
    
17.         c(i)=-1-(h/2)*p(x);8 h\" }& |% T; T- H' J* h
    
18.         d(i)=-h*h*r(x);
    - F$ O3 p, I  S. }
19.      end
    : ?# F+ x3 [  a& y. D2 ~
20.      x=b0-h;/ K- Z  Y4 {$ t( \4 z% @# z
    
21.      a(N)=2+h*h*q(x);1 y, c* h1 m% B! h- N& v& T* u
    
22.      c(N)=-1-(h/2)*p(x);1 a5 B! S& T7 E: Z! I8 q
    
23.      d(N)=-h*h*r(x)+(1-(h/2)*p(x))*bt;4 P) @) ]3 W- s/ i& x5 R0 x+ M* I
    
24.      %%%%%%%%%追赶法%%%%%%%%%%%%%%%%%%
    ' z, w% L# z& z$ A6 m  V; l( n
25.      %y=trisys(c,a,b,d)
    6 Y- d0 Q1 f  F# ^. O
26.      L(1)=a(1);
    4 ^3 q; B+ b. _
27.      u(1)=b(1)/a(1);4 H\" t' G* {! E' f! A8 H+ W7 U
    
28.      for i=2:N-17 a  f% f7 d: m2 p0 a
    
29.          L(i)=a(i)-c(i)*u(i-1);8 a, J9 k3 _$ I* E
    
30.          u(i)=b(i)/L(i);
    5 N% f% z4 b  p  A2 C
31.      end
    , p) E\" I2 h2 B4 E
32.      L(N)=a(N)-c(N)*u(N-1);
    & k5 n\" b$ p7 Q' T3 N8 G6 u
33.      z(1)=d(1)/L(1);# r4 q, Z( V0 ~1 @
    
34.      for i=2:N4 c8 Y2 i1 ]% n, m. F6 q* d1 Y
    
35.          z(i)=(d(i)-c(i)*z(i-1))/L(i);9 g# \4 B% u  v: V\" m& T
    
36.      end
    + ?% o\" S8 E7 \& M
37.      y(N)=z(N);
    . i% [, O7 ~/ r) I
38.      for i=N-1:-1:14 A! m4 V# c% V; q& s* O; g- k
    
39.          y(i)=z(i)-u(i)*y(i+1);
    / w4 j9 e2 M7 S- b7 k
40.      end1 ]& \* |, m( e' ]# U
    
41.      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    $ K, s1 {# {! W) A. C5 |9 b
42.      Y=[af,y,bt];
    9 M5 {/ y  \. Y0 M; C\" ?- F' V! b
43. for i=1:N+2( Z+ [\" p$ Y% s8 Q+ z3 {
    
44.       x=a0+(i-1)*h;( _7 _1 U, Q* }  |- H6 k& \& @
    
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;3 z7 P( m) p1 i
    
46. end
    # D. [, I; ]$ l8 V) \4 H* Z6 E: j
47. disp('下面两列分别是数值解和近似解');
    7 A1 J- q) f\" I' j0 X
48. re=[Y'      zj']

复制代码

/ ^4 y8 p0 Q: c