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

2744557306

使用有限差分法和托马斯算法（或追赶法）对一个二阶线性边值问题进行数值求解。这种方法通常用于数值解微分方程。
以下是代码的简要解释：
0 }; i% N1 E9 P$ E0 [
1.使用 inline 函数定义了三个函数 p(x)、q(x) 和 r(x)，它们表示微分方程的系数。
2.设置了参数，如间隔数 N、初始和边界条件 a0、b0、af、bt 以及间隔大小 h。
, V+ _3 G: p# M( `% N3 G& [3.基于微分方程的有限差分离散化，计算了系数 a、b、c 和 d。
; y  h9 a) O# [4.使用托马斯算法（或追赶法）解决了三对角方程组。
5.将结果与由数组 zj 表示的解析解进行了比较。
6.将数值解和解析解并排显示，以便比较。

1. p=inline('-2/x');
   4 g4 q$ o% f5 a: E9 W
2.   q=inline('2/x^2');
   $ m\" W. F% H# u  s+ v
3.   r=inline('sin(log10(x)/log10(exp(1)))/x^2');) r. U  N' J' a5 J( K2 w
   
4.   N=9;
   2 ]+ @2 o9 q! N0 t
5.   a0=1;b0=2;
   & [1 ~' p' e7 n' Q' H
6.   af=1;bt=2;' m5 f3 N+ q, ?
   
7.   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%* X% U/ S# x1 m0 j
   
8.   h=(b0-a0)/(N+1);
   . S* B6 _* z\" q' Q\" j
9.   x=a0+h;& W8 R* P$ ?* Q, Z
   
10.   a(1)=2+h*h*q(x);
    , K% g6 U& l7 \* Q$ ]4 z\" O6 l: {
11.   b(1)=-1+(h/2)*p(x);& S; O0 I8 ^. t, v) _; e
    
12.   d(1)=-h*h*r(x)+(1+(h/2)*p(x))*af;
    7 g( N; z\" U6 k
13.     for i=2:N-1  S1 K- X+ i( N, x1 A
    
14.         x=a0+i*h;
    ! [9 H; F1 Y# G, E: N
15.         a(i)=2+h*h*q(x);
    . u& G3 R$ Q# p: ?\" d& M6 L- s
16.         b(i)=-1+(h/2)*p(x);
    9 A! A, K0 e$ U1 F7 l2 v
17.         c(i)=-1-(h/2)*p(x);
    1 u; j& T: T0 d& L  B: c
18.         d(i)=-h*h*r(x);
    3 N! t1 n; J/ ~4 L8 Y
19.      end
    * N) r* D& j6 F8 f( E+ v; e+ S
20.      x=b0-h;/ H# k# {/ F( }# _2 ~* T- ?* Q
    
21.      a(N)=2+h*h*q(x);; ]! o* m( f9 a' F; Q  j
    
22.      c(N)=-1-(h/2)*p(x);
      H7 \+ D/ V5 w
23.      d(N)=-h*h*r(x)+(1-(h/2)*p(x))*bt;
    ! q# C\" H8 X8 n) S/ f5 I
24.      %%%%%%%%%追赶法%%%%%%%%%%%%%%%%%%
    ( _0 K# O8 t# I, l& x
25.      %y=trisys(c,a,b,d)- W7 U' y# e! D# Y5 ]$ o. B5 e' ^
    
26.      L(1)=a(1);
    4 b* N. p2 S: P9 O9 s
27.      u(1)=b(1)/a(1);
    7 L5 t\" s: O& ?/ n\" u4 S
28.      for i=2:N-19 F\" t8 R8 C) r: k2 b# q+ O, R\" A8 l
    
29.          L(i)=a(i)-c(i)*u(i-1);' f$ {7 _, R  q) b: n
    
30.          u(i)=b(i)/L(i);7 O5 `) d) l  F) ]
    
31.      end( ?1 m, j* `. L
    
32.      L(N)=a(N)-c(N)*u(N-1);2 y9 Q& M- C$ k8 O. A& Q
    
33.      z(1)=d(1)/L(1);
    / o, p- s1 j3 M9 f: @+ a/ t: ^
34.      for i=2:N
    * o4 c$ U) @1 b\" \$ W$ C
35.          z(i)=(d(i)-c(i)*z(i-1))/L(i);, g% F% t% J. a7 T$ {. y6 }4 N\" v
    
36.      end4 }; w1 d: J1 B' \, y
    
37.      y(N)=z(N);5 h  c7 E\" v- q6 s# ~0 @
    
38.      for i=N-1:-1:1
    4 h2 X0 s( \& \\" h; m5 v7 o/ K
39.          y(i)=z(i)-u(i)*y(i+1);4 h8 n. L; O5 F2 \. V  b
    
40.      end& R7 w2 f' i; w  h
    
41.      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%: k  w+ v$ x$ j! Q7 k) c
    
42.      Y=[af,y,bt];
    $ _* k/ m; l3 D9 f  p+ S
43. for i=1:N+2; Y2 ^: i: `) }% B. T* L
    
44.       x=a0+(i-1)*h;/ [% O9 ~, C, P5 m5 e\" U- y. Q
    
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;
    , `2 P$ p9 W% M$ H6 p* W( Q& s% [
46. end
    ! {! ]6 g% W; X. R2 D4 k+ _
47. disp('下面两列分别是数值解和近似解');
    4 @' `( n& M1 a1 C
48. re=[Y'      zj']

复制代码