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

[复制链接]

字体大小: 正常放大

2744557306

1192 主题	4 听众	2946 积分

该用户从未签到

电梯直达

1^#

发表于 2024-1-3 09:34 |只看该作者 |倒序浏览

|招呼Ta 关注Ta

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

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

p=inline('-2/x');: z/ B\" m! r' {' e: ^$ I
q=inline('2/x^2');1 |8 K; C5 i3 d* i( n6 R ^3 A4 k8 o) o
r=inline('sin(log10(x)/log10(exp(1)))/x^2');
6 [( s E% F# I& K4 h8 o/ s
N=9;; B8 _4 ^' P, I. l8 I; G
a0=1;b0=2;3 W4 U$ v7 x% ^
af=1;bt=2;% i9 I8 M' P( W) ?- |- {4 D' H
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
! o( n( I& C\" h% Y\" s2 |- `/ c
h=(b0-a0)/(N+1);
5 _0 Q0 ?6 ]/ P; q3 Y( S
x=a0+h;- D* c8 K\" J\" l7 A
a(1)=2+h*h*q(x);$ }9 n0 B- U6 M) ~: ?' R3 p
b(1)=-1+(h/2)*p(x);7 }/ |$ i$ Z* h. C r\" O* N% H
d(1)=-h*h*r(x)+(1+(h/2)*p(x))*af;
# `8 U* N- `/ [' v. n j/ N; o
for i=2:N-1
6 }/ h+ z* H# f; l/ y# J
x=a0+i*h;0 `. w\" m2 |# A t+ k
a(i)=2+h*h*q(x);
/ B5 [0 r3 b, h7 _& g. P& Q
b(i)=-1+(h/2)*p(x);* \4 x3 z$ d2 ~2 C9 h. r: E
c(i)=-1-(h/2)*p(x);/ }# d& O$ {. O$ }) g( q( u/ o
d(i)=-h*h*r(x);! M6 J' c; Y) q+ \\" S4 M
end+ t {: A9 j. Y$ @\" r: c
x=b0-h;. W3 }3 l- E0 A# {
a(N)=2+h*h*q(x);0 J! A- H* b/ K* b1 O; m
c(N)=-1-(h/2)*p(x);3 P! O$ a% d+ f- n4 O
d(N)=-h*h*r(x)+(1-(h/2)*p(x))*bt;5 s! f' ?- @! i2 P- h
%%%%%%%%%追赶法%%%%%%%%%%%%%%%%%%9 o3 e+ a/ ^. C
%y=trisys(c,a,b,d)
- Q$ I% Q' ~6 N4 c5 K
L(1)=a(1);
+ j, a! V8 B2 {' ]! N$ V
u(1)=b(1)/a(1);7 t- R! O* h3 H1 b
for i=2:N-15 H! J+ _7 g6 I
L(i)=a(i)-c(i)*u(i-1);$ r( g$ n/ X; q+ ]3 m( B9 _
u(i)=b(i)/L(i);
& ^0 l5 h. a, l& M7 a
end7 b$ V8 H. p3 M8 \7 J7 p( W1 z
L(N)=a(N)-c(N)*u(N-1);# [, ~7 j, f4 d
z(1)=d(1)/L(1);
0 G3 L+ d& K1 ` t7 n6 ?
for i=2:N
4 y) M+ Y# j8 p$ t0 A, [
z(i)=(d(i)-c(i)*z(i-1))/L(i);4 Z3 U8 c6 L9 |8 J
end
* O( w) m6 _5 B* g
y(N)=z(N);
/ w' P5 ~4 s. A
for i=N-1:-1:1, `; C* a) z* O
y(i)=z(i)-u(i)*y(i+1);* F- R# L; q F7 W
end
1 G& D: X$ [3 R\" r+ Y
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
( _: [7 N8 ^, V( E' o: W
Y=[af,y,bt];
\: G\" S& X+ t# V# a: h
for i=1:N+2
2 Q% m( {/ U, ~& t9 ^% h- j\" m
x=a0+(i-1)*h;
3 g# E6 D/ x2 \% E8 h+ t\" O
zj(i)=1.1392070132*x-0.03920701320/x^2-3*sin(log10(x)/log10(exp(1)))/10-cos(log10(x)/log10(exp(1)))/10;
6 O1 a9 k' b& k( a6 p& g$ s1 @
end
$ H: p) ^3 d: p7 n
disp('下面两列分别是数值解和近似解');
# X( j; M% c. {
re=[Y' zj']