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

2744557306

使用有限差分法和托马斯算法（或追赶法）对一个二阶线性边值问题进行数值求解。这种方法通常用于数值解微分方程。
! [4 |: Z% C' K# [6 `以下是代码的简要解释：
! F$ _9 o  v0 I5 P: ?6 Z
1.使用 inline 函数定义了三个函数 p(x)、q(x) 和 r(x)，它们表示微分方程的系数。
$ k0 V' B$ J8 H0 C' z2.设置了参数，如间隔数 N、初始和边界条件 a0、b0、af、bt 以及间隔大小 h。
: y7 _7 D! t" O" H; `8 V+ x. S3.基于微分方程的有限差分离散化，计算了系数 a、b、c 和 d。
4.使用托马斯算法（或追赶法）解决了三对角方程组。
5.将结果与由数组 zj 表示的解析解进行了比较。
& d0 T" f& P" W) ]% O+ `6 S1 N( F, @6.将数值解和解析解并排显示，以便比较。

1. p=inline('-2/x');8 R% P2 K+ M4 y9 d, ]5 i/ G# Q
   
2.   q=inline('2/x^2');+ [/ P9 y( ~. a$ N& {; B# [
   
3.   r=inline('sin(log10(x)/log10(exp(1)))/x^2');3 Q: ?& p! s( V' \. `
   
4.   N=9;6 d( d\" L$ }) `  ]6 Z
   
5.   a0=1;b0=2;
   + I( v3 G$ ^* O% N, @
6.   af=1;bt=2;& q8 {& N9 i) `' ~
   
7.   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   4 h1 p2 @& x; O* i2 O( \
8.   h=(b0-a0)/(N+1);
   # R2 \$ w9 A* \2 h% p  o8 W
9.   x=a0+h;+ H. n. p6 G; A+ g& ?
   
10.   a(1)=2+h*h*q(x);- U' I) d  {$ U' a! N/ F- e# w% J
    
11.   b(1)=-1+(h/2)*p(x);
    1 {* m- w+ S/ L5 u\" l/ X! X9 \
12.   d(1)=-h*h*r(x)+(1+(h/2)*p(x))*af;
    1 \3 c; I& P: \; K  [9 Z! X
13.     for i=2:N-1- p! e! q6 \5 k6 O2 I6 s! s
    
14.         x=a0+i*h;
    9 [6 g7 E% X0 v, N' \
15.         a(i)=2+h*h*q(x);
    / \1 D: m$ n! Y) x/ L8 y1 J
16.         b(i)=-1+(h/2)*p(x);
    \" v0 |5 w; z: O5 f/ O3 @
17.         c(i)=-1-(h/2)*p(x);
    0 Q( r( R/ V7 E1 M  b; R
18.         d(i)=-h*h*r(x);6 t' S# {# A; R) ~, _3 U
    
19.      end
    \" a% s9 c8 `\" \2 N
20.      x=b0-h;1 Q2 S$ Q7 t\" m1 ^8 m/ W0 G
    
21.      a(N)=2+h*h*q(x);
    & R8 }) t$ G\" C+ Y: s
22.      c(N)=-1-(h/2)*p(x);4 u; c2 _* {, F7 H
    
23.      d(N)=-h*h*r(x)+(1-(h/2)*p(x))*bt;
    2 ?2 |* A; N1 g# s6 T9 V- J# @$ X
24.      %%%%%%%%%追赶法%%%%%%%%%%%%%%%%%%5 }& k\" u5 o6 {7 {3 m$ m8 _2 O
    
25.      %y=trisys(c,a,b,d)
    2 i! U/ m' h# h, H: A  }
26.      L(1)=a(1);
    ( s% d$ r& ^\" w8 Q
27.      u(1)=b(1)/a(1);3 l) Y# T; s/ q! B9 B
    
28.      for i=2:N-1
    / E8 e, V; R3 O) t9 Y
29.          L(i)=a(i)-c(i)*u(i-1);
    ; s3 u+ ^$ q7 `' Y% R( j
30.          u(i)=b(i)/L(i);
    % f0 _4 \0 s4 g& m, c2 C) G
31.      end
    ; H+ y$ i$ G( L- U' J: v
32.      L(N)=a(N)-c(N)*u(N-1);2 ^6 p0 [2 k- N5 \
    
33.      z(1)=d(1)/L(1);8 m2 K  \5 `1 x
    
34.      for i=2:N
    / N9 |! V, i  e# P5 M
35.          z(i)=(d(i)-c(i)*z(i-1))/L(i);
      e  y  U, t\" r8 D
36.      end
    / e5 F8 L0 z% Z# k
37.      y(N)=z(N);
    # [1 L$ s( }% v1 H# t: R
38.      for i=N-1:-1:1
    0 \9 S/ i$ [  h0 V5 S+ d
39.          y(i)=z(i)-u(i)*y(i+1);
    $ l6 K7 d0 y( b) B
40.      end; L9 C, ?7 J( _
    
41.      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    * D7 `- b2 t. U# _; u
42.      Y=[af,y,bt];
    7 N2 ^\" q4 U( H# X
43. for i=1:N+2
    7 j# o4 V, R$ J8 Z' }% ]+ g
44.       x=a0+(i-1)*h;
      l$ |1 A* X- O- i$ B2 E+ n; E) j
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;: |% m/ {$ \- w* u. X7 A+ ]& Z4 J1 w
    
46. end+ u# o. E  Q5 O& `
    
47. disp('下面两列分别是数值解和近似解');
    3 O  r; `\" L3 }2 d
48. re=[Y'      zj']

复制代码