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

[复制链接]

字体大小: 正常放大

2744557306

1175 主题	4 听众	2818 积分

该用户从未签到

电梯直达

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');
/ M6 O6 W$ {! I6 t8 t& O\" d
q=inline('2/x^2');
% A0 b- S; O4 [' \
r=inline('sin(log10(x)/log10(exp(1)))/x^2');2 T4 Y1 @+ F. H; n8 C( h+ J
N=9;* T. [7 r3 f* W+ s; h# \
a0=1;b0=2;
- h/ J& g+ H- D% R6 o2 A
af=1;bt=2;9 C n% U5 ^0 Y% [% J, W
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
. Q) V4 y! p/ w
h=(b0-a0)/(N+1);
& N4 T8 X: ~2 [\" L
x=a0+h;3 M( y( C. v% P! U1 ~& M3 z
a(1)=2+h*h*q(x);
8 u. Y7 f' A: a' s' g1 {
b(1)=-1+(h/2)*p(x); C; _8 a, T, X' `. ~$ `
d(1)=-h*h*r(x)+(1+(h/2)*p(x))*af;, X6 V/ I+ a: e
for i=2:N-11 c0 O8 L) p* L; b: T4 A e3 i3 _1 S
x=a0+i*h;
* b+ j1 x6 v8 A# p# z0 ~
a(i)=2+h*h*q(x);
+ v& b3 G1 o) \* F& @
b(i)=-1+(h/2)*p(x);
f* P/ o: t; b: S
c(i)=-1-(h/2)*p(x);
# ], b- \ H% a- V: q* x4 ?+ R9 C
d(i)=-h*h*r(x);( r4 g7 F g. C) n
end
) `. H8 [8 ~5 K* U0 ~- V
x=b0-h;
4 V: `9 A c\" N7 a' \
a(N)=2+h*h*q(x);0 E8 @! H+ u& x\" c; l/ V
c(N)=-1-(h/2)*p(x);
( M$ F% c5 @) s: ?3 e3 N
d(N)=-h*h*r(x)+(1-(h/2)*p(x))*bt;: P T\" W( _* `. C h9 D. n
%%%%%%%%%追赶法%%%%%%%%%%%%%%%%%%$ \( ]$ @. h9 ~
%y=trisys(c,a,b,d)8 K+ N, u% J# e. l x* p0 V
L(1)=a(1);4 J; a3 b1 z' } }+ ~
u(1)=b(1)/a(1);2 n( K1 k3 D) L V
for i=2:N-11 ~9 X\" a) w$ E+ z1 b/ w* [
L(i)=a(i)-c(i)*u(i-1);
: t) b; }0 h7 g) p4 P/ N
u(i)=b(i)/L(i);
p0 {, {* X$ S& S- G) E
end# p( s3 v% y$ g
L(N)=a(N)-c(N)*u(N-1);
: v; W8 C) M9 h: m\" @: C
z(1)=d(1)/L(1);9 ~4 x+ G, g- ]$ w
for i=2:N
% p# |/ \# N. G/ ~, g* T
z(i)=(d(i)-c(i)*z(i-1))/L(i);/ j) Y4 U$ t: w# c7 C
end) {( b, Z( N) i\" v8 E/ e9 s: \
y(N)=z(N);
* m7 m2 \7 x% H\" n: u; U
for i=N-1:-1:1) A9 d- r* |: H/ |- @\" o
y(i)=z(i)-u(i)*y(i+1);
5 H! g\" z8 k* `% J# G# I! G% b f
end. R8 `. |4 Z' Q' r. u
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%6 J) h: W+ w3 ^$ {3 T+ ]
Y=[af,y,bt];\" v, D7 @( `7 w4 ]% X
for i=1:N+2
6 {# o' P, O F' p0 E
x=a0+(i-1)*h;
2 j; l5 }6 j) d4 w) c6 d
zj(i)=1.1392070132*x-0.03920701320/x^2-3*sin(log10(x)/log10(exp(1)))/10-cos(log10(x)/log10(exp(1)))/10;
+ J. P2 {& i8 y8 V. N6 Q7 X0 l
end0 i- B0 |, h+ V
disp('下面两列分别是数值解和近似解'); }+ w; @$ |% l W3 K/ V$ _
re=[Y' zj']