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

2744557306

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

1.使用 inline 函数定义了三个函数 p(x)、q(x) 和 r(x)，它们表示微分方程的系数。
; v7 X$ p5 j: r2.设置了参数，如间隔数 N、初始和边界条件 a0、b0、af、bt 以及间隔大小 h。
3.基于微分方程的有限差分离散化，计算了系数 a、b、c 和 d。
4.使用托马斯算法（或追赶法）解决了三对角方程组。
5.将结果与由数组 zj 表示的解析解进行了比较。
0 K( y# V, ]' k& O3 Y; \6.将数值解和解析解并排显示，以便比较。

1. p=inline('-2/x');2 h& w7 P# a) K5 _, B
   
2.   q=inline('2/x^2');! ], E' R+ o4 j5 C3 h
   
3.   r=inline('sin(log10(x)/log10(exp(1)))/x^2');
   ( M# D  }' i' y2 [
4.   N=9;
   , L\" y; Q% p' Q, ^
5.   a0=1;b0=2;
   \" T( ]6 i, s9 g, O2 K! O6 T
6.   af=1;bt=2;
   $ a9 h( B: c6 K8 y6 ^
7.   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   7 V# `$ G\" `% r- X: f
8.   h=(b0-a0)/(N+1);
   \" l3 |; e0 y# z\" ^  r  Y3 d5 |
9.   x=a0+h;
   7 v/ {; p% E, Q8 I4 g9 h# s+ m6 O7 V/ f
10.   a(1)=2+h*h*q(x);# m# K2 `7 q% r; L2 _& G3 c' J. C# b0 q
    
11.   b(1)=-1+(h/2)*p(x);9 c9 f/ G1 J0 x7 l& R# W9 ?* d& y
    
12.   d(1)=-h*h*r(x)+(1+(h/2)*p(x))*af;3 j( e9 Z$ w& ], d2 F7 p
    
13.     for i=2:N-1
    ; l3 ^4 ]# E0 K5 c. V
14.         x=a0+i*h;# d( k5 h8 I0 o4 z; |9 c7 t# j
    
15.         a(i)=2+h*h*q(x);+ @: O- M$ r\" d& M- S$ ]
    
16.         b(i)=-1+(h/2)*p(x);
    0 j( T, @, ], ]! {/ @: ]
17.         c(i)=-1-(h/2)*p(x);
    6 e$ ]6 ^3 R5 a  j
18.         d(i)=-h*h*r(x);
    - A) H6 I8 d+ r3 h0 Y1 x
19.      end1 m/ d/ o5 G\" U: S/ G
    
20.      x=b0-h;
    ; C5 w7 y3 U# |+ z  Y* V
21.      a(N)=2+h*h*q(x);4 g: V- p5 ?5 {: u( K
    
22.      c(N)=-1-(h/2)*p(x);
      l; n8 N. I  t: @\" J/ Q5 z
23.      d(N)=-h*h*r(x)+(1-(h/2)*p(x))*bt;5 u% G. \! E' J' j9 D) C$ C
    
24.      %%%%%%%%%追赶法%%%%%%%%%%%%%%%%%%4 I5 X# m% p/ q. J2 O+ L$ O; I3 i
    
25.      %y=trisys(c,a,b,d)' W% A2 k% D$ H
    
26.      L(1)=a(1);
    $ M6 d2 L* T\" }
27.      u(1)=b(1)/a(1);
    ) K/ ~( Q8 }! k\" j
28.      for i=2:N-1
    . L+ `: k& w* j4 t9 ^' W$ U
29.          L(i)=a(i)-c(i)*u(i-1);
    ) L* I! J! b! J8 t
30.          u(i)=b(i)/L(i);/ A; W$ n6 W  V& e; F+ \
    
31.      end
    - i: v% }% f; P2 ]& L0 ~2 U\" B) ]
32.      L(N)=a(N)-c(N)*u(N-1);! p! E0 A  G\" d$ L: a
    
33.      z(1)=d(1)/L(1);
    5 [, `2 d, {6 b8 U; h
34.      for i=2:N
    , w0 l/ g; c9 N$ j1 M, }
35.          z(i)=(d(i)-c(i)*z(i-1))/L(i);* n1 @: r3 A  y; H9 W
    
36.      end) M0 K0 ~) F' w1 m\" F8 J7 B4 u
    
37.      y(N)=z(N);$ e  i# {\" B% Z9 W( y
    
38.      for i=N-1:-1:1
    3 D9 [& i- o% H4 W
39.          y(i)=z(i)-u(i)*y(i+1);
    & W: `8 [6 B& c% q7 b\" V
40.      end
    5 L9 B4 e2 X! `% D\" \9 T6 O
41.      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    ; l2 C- X) P! m8 R; X, |
42.      Y=[af,y,bt];
    5 Q$ L: X$ w8 T' h+ ^
43. for i=1:N+2' \. p% [! D4 w* s5 a
    
44.       x=a0+(i-1)*h;9 m7 U0 ]  a- H) q* h5 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;
    , n$ X( }# d9 H\" k! p8 U/ c
46. end
      X9 X1 N6 n( E7 b# A. ]- h\" E' _0 A\" n
47. disp('下面两列分别是数值解和近似解');/ K# x8 a2 z$ L( I& l  p
    
48. re=[Y'      zj']

复制代码