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

2744557306

使用有限差分法和托马斯算法（或追赶法）对一个二阶线性边值问题进行数值求解。这种方法通常用于数值解微分方程。
以下是代码的简要解释：
' @) w4 \% t! N
1.使用 inline 函数定义了三个函数 p(x)、q(x) 和 r(x)，它们表示微分方程的系数。
2.设置了参数，如间隔数 N、初始和边界条件 a0、b0、af、bt 以及间隔大小 h。
3 x2 j# P( M; ~/ \; c0 m3.基于微分方程的有限差分离散化，计算了系数 a、b、c 和 d。
( ~) K6 h) m( D, J- T4.使用托马斯算法（或追赶法）解决了三对角方程组。
5.将结果与由数组 zj 表示的解析解进行了比较。
2 f& x, Y1 `( E( A  A; o6.将数值解和解析解并排显示，以便比较。

1. p=inline('-2/x');
     [) t3 s- \' u/ e2 n\" e3 W
2.   q=inline('2/x^2');
   9 Q% B+ A# ?9 r6 a( W# }  H4 o8 L- c
3.   r=inline('sin(log10(x)/log10(exp(1)))/x^2');
   * D& X% Y' J1 Q' j( U' l* f
4.   N=9;: F! F# o  `$ X  F: c! R
   
5.   a0=1;b0=2;
   ! G7 a1 B\" `0 C
6.   af=1;bt=2;
   9 A# p  h8 Y& `$ X6 U+ x! q
7.   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5 q& L  V) `8 A' m3 F' v
   
8.   h=(b0-a0)/(N+1);
   # ]  h! X  b1 b8 ]6 H
9.   x=a0+h;
   5 Y, W  \* i- W
10.   a(1)=2+h*h*q(x);
    ; L\" i9 A' [- H5 e- V8 O5 i  ^
11.   b(1)=-1+(h/2)*p(x);
    : Q3 C) P2 a; ~
12.   d(1)=-h*h*r(x)+(1+(h/2)*p(x))*af;( N4 q\" g2 p( P0 r. J
    
13.     for i=2:N-1- S: w6 _+ g\" U5 C6 g$ X
    
14.         x=a0+i*h;
    3 Q+ U8 d- s' w5 [+ x' q9 r6 u
15.         a(i)=2+h*h*q(x);
    ! H6 l. t  F6 g! [3 u2 v  z( `
16.         b(i)=-1+(h/2)*p(x);
    , B( X- g) g% ]- O\" F: i
17.         c(i)=-1-(h/2)*p(x);
    ! J3 ?# _. ]$ j. v) \% R$ m
18.         d(i)=-h*h*r(x);
    6 Q+ S3 X! y/ l2 z1 w2 K
19.      end1 E( Z  |2 }1 I1 b- W5 _( b+ r
    
20.      x=b0-h;
    , N5 g7 v! T3 X0 z
21.      a(N)=2+h*h*q(x);
    . N# q2 ^9 q6 j' }! t3 K. x
22.      c(N)=-1-(h/2)*p(x);* Y2 K+ L# j* c+ x& [
    
23.      d(N)=-h*h*r(x)+(1-(h/2)*p(x))*bt;
    , ?) Z$ j3 @+ g& |0 y# m, ~4 I+ z
24.      %%%%%%%%%追赶法%%%%%%%%%%%%%%%%%%
    : y, v: |$ v7 {3 |8 y; m
25.      %y=trisys(c,a,b,d)
    & y: ^/ h+ p4 H+ i, j& [3 K
26.      L(1)=a(1);) i) G+ y\" ^* m9 i
    
27.      u(1)=b(1)/a(1);
    : C. Y/ M4 E# @
28.      for i=2:N-1
    5 U  c- a& S/ K
29.          L(i)=a(i)-c(i)*u(i-1);
    2 A) Z- s# s\" I0 v) d  _7 s; Z
30.          u(i)=b(i)/L(i);
    2 }$ ?8 R' t5 Q
31.      end
    / }\" K( @/ ]* G+ i4 O
32.      L(N)=a(N)-c(N)*u(N-1);
    5 X$ \( e8 m/ g# _
33.      z(1)=d(1)/L(1);, s6 a! {# `# ~$ N
    
34.      for i=2:N
    6 z! n# Z5 ^( ]! c. D
35.          z(i)=(d(i)-c(i)*z(i-1))/L(i);
    % {* d6 ]4 Z* J% b% q. M- r
36.      end! w3 ]6 I5 U1 N* ^5 M
    
37.      y(N)=z(N);
    2 f/ d# z2 E6 L+ K4 y
38.      for i=N-1:-1:14 T9 ^' S# @' ~7 @# q
    
39.          y(i)=z(i)-u(i)*y(i+1);* ~5 V# o4 ?: h6 m; Z% ]1 |4 w
    
40.      end
    9 q0 ~* b+ |; Z
41.      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    5 d8 P- _/ E  ^
42.      Y=[af,y,bt];$ V' Q/ c1 H, B
    
43. for i=1:N+2\" t/ M+ k7 z% c  q$ i, j' ^
    
44.       x=a0+(i-1)*h;9 G' ]7 E% O2 Y. E9 C, v8 {
    
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;
    ; }1 ?; M; A' k* C
46. end
    7 G9 B2 U4 V! o0 t8 O
47. disp('下面两列分别是数值解和近似解');! B' d  F% K7 c, `4 d1 x
    
48. re=[Y'      zj']

复制代码

9 `3 y4 W$ \2 `