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

[复制链接]

字体大小: 正常放大

2744557306

1186 主题	4 听众	2922 积分

该用户从未签到

电梯直达

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');+ B& |/ g6 {' j3 y
q=inline('2/x^2');
5 \2 X3 F4 A8 B$ F/ v
r=inline('sin(log10(x)/log10(exp(1)))/x^2');! }9 x8 M$ z2 f# W7 {
N=9;
4 ]5 `) t% ?) F- m6 F\" ?5 n' o
a0=1;b0=2;+ K4 f6 X6 v; }6 j K
af=1;bt=2;
5 k; w: t; D2 Z( D: x, P
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
* U6 d9 ?+ ~2 t6 l0 Q
h=(b0-a0)/(N+1);# O( C# M1 w% l# p5 X3 N( N+ Z# C
x=a0+h;/ p% `1 o8 x# u6 h1 z
a(1)=2+h*h*q(x);$ T! Y6 J b0 t- `/ P
b(1)=-1+(h/2)*p(x);
. _. r k# p2 {
d(1)=-h*h*r(x)+(1+(h/2)*p(x))*af;9 ^* b/ x0 _, P' ?3 h
for i=2:N-1
$ U# G1 @' O- P* N3 M$ P
x=a0+i*h;
6 X5 `9 o/ H: H4 g4 |$ R$ |1 W\" N
a(i)=2+h*h*q(x);. [* U! j\" F, y5 _
b(i)=-1+(h/2)*p(x);9 I+ \7 g7 _7 y0 ] h4 B1 G; F
c(i)=-1-(h/2)*p(x);: Q: s( B! `; J; I
d(i)=-h*h*r(x);
) b! D$ Y+ O* F! v8 `
end d8 @& \. M \* ?# k/ e
x=b0-h;& H, j3 T/ m i. F9 e! t- Q
a(N)=2+h*h*q(x);4 e# L1 _, \; Z2 n
c(N)=-1-(h/2)*p(x);# a& F- O; o4 W9 s+ l0 o6 x
d(N)=-h*h*r(x)+(1-(h/2)*p(x))*bt;
7 v( ^' U/ L0 W
%%%%%%%%%追赶法%%%%%%%%%%%%%%%%%%0 h: p1 V7 O) w; J, ~: W
%y=trisys(c,a,b,d)9 X9 H0 p- z/ n9 p+ O( p
L(1)=a(1);1 y+ v# M\" t3 W5 B; p$ I\" o& e* L
u(1)=b(1)/a(1);
$ l/ |$ T6 y$ ^+ [: W; ~
for i=2:N-14 ?3 r; l\" V: W9 E
L(i)=a(i)-c(i)*u(i-1);
`, {/ u% k& r v9 M
u(i)=b(i)/L(i);
4 Q% u: h% @' }3 G
end' D& v0 U! J% U% H1 O! L
L(N)=a(N)-c(N)*u(N-1);! W' n& {! R; L2 Z ], R
z(1)=d(1)/L(1);5 O8 p5 S2 f/ ?& Y |9 t) O6 n( M. {
for i=2:N) \) \7 A1 M y
z(i)=(d(i)-c(i)*z(i-1))/L(i);
/ N\" E7 y' |+ f
end4 k+ f+ C- K, ^, E1 Y% R y$ R
y(N)=z(N);
8 F- L& `/ e6 ^3 n1 d5 L' u8 A
for i=N-1:-1:16 d: S4 I& g; \/ j
y(i)=z(i)-u(i)*y(i+1);
+ T& d) Y\" Q% ]( b
end
+ l% N4 `3 M: I) Q0 o! N& N* P' T
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%$ v2 t5 M2 o. Y7 N
Y=[af,y,bt];4 ~1 G. c4 p: \' N
for i=1:N+2
/ y6 K% ~' g, M
x=a0+(i-1)*h;
9 j! e H2 f) W\" f
zj(i)=1.1392070132*x-0.03920701320/x^2-3*sin(log10(x)/log10(exp(1)))/10-cos(log10(x)/log10(exp(1)))/10;/ R. a6 P. C, }* b
end
5 m6 `) _ a) m, }/ ^
disp('下面两列分别是数值解和近似解');
( F\" y/ P' C. s( c. W
re=[Y' zj']