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

[复制链接]

字体大小: 正常放大

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');+ L\" A/ V8 A3 d+ e: Q0 r: a: }4 r
q=inline('2/x^2');
& b* C; ?% Y* _7 B, v, W
r=inline('sin(log10(x)/log10(exp(1)))/x^2');
& @+ u! |; n N
N=9;- p8 g\" }- u' x( ]; ^$ {; _; L
a0=1;b0=2;) _0 }\" y* u7 D' H1 d
af=1;bt=2;9 M/ ^% d8 T g/ }* c- O
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
. u5 {9 T1 _0 m( E+ w
h=(b0-a0)/(N+1);. |0 q5 S0 j: J4 ]* Q! F) Q
x=a0+h;
1 x t! I/ d' c. A
a(1)=2+h*h*q(x);
8 c# I/ ?: q, G; V# I
b(1)=-1+(h/2)*p(x);
7 N9 Z& B4 z. s/ ?( h
d(1)=-h*h*r(x)+(1+(h/2)*p(x))*af;2 L, |: r1 G/ L$ Z9 m
for i=2:N-1\" \# ]# ]* N* F& i* o- V
x=a0+i*h;4 h) L\" z2 ?\" l2 _8 a1 z
a(i)=2+h*h*q(x);
N9 F9 M9 c) H\" e& A4 d
b(i)=-1+(h/2)*p(x);( ~9 D$ n$ b\" L7 ]
c(i)=-1-(h/2)*p(x);5 m4 L; C! y8 w5 [; D( \2 v\" A
d(i)=-h*h*r(x);
& K6 D, i& _/ _\" L2 S) v
end
* ?- n7 g6 z# h/ ^. y% s- R
x=b0-h;
7 o1 {# z- y' }* ^6 h
a(N)=2+h*h*q(x);5 k' h P9 n5 N( \4 ~
c(N)=-1-(h/2)*p(x);/ B. a& ]\" S2 M/ f0 l7 r! Q
d(N)=-h*h*r(x)+(1-(h/2)*p(x))*bt;
1 E2 R# q1 B) {4 r
%%%%%%%%%追赶法%%%%%%%%%%%%%%%%%%
! g& y9 j4 S9 L! C\" q# U! W9 o' y
%y=trisys(c,a,b,d)
5 ?/ O, \! S+ z& u* y3 ~
L(1)=a(1);, s: @\" B3 T. _! r' `8 @
u(1)=b(1)/a(1);0 V6 k, i( `3 o8 l$ Z6 t
for i=2:N-1/ ~' ^/ D* H' T, ` m. C
L(i)=a(i)-c(i)*u(i-1);5 j- B: N( H3 w9 o! C
u(i)=b(i)/L(i);
0 y$ @# D2 Z& h4 z
end
4 w' ~( K6 B* v3 T' Q
L(N)=a(N)-c(N)*u(N-1);$ ], |' h3 Y0 ]' R5 {0 U1 g
z(1)=d(1)/L(1);9 s! }; @7 S% H
for i=2:N
# x8 o% o9 b& ]8 s/ _+ V% g
z(i)=(d(i)-c(i)*z(i-1))/L(i);6 ^$ d. @ @& y5 H! k2 Z
end
- o/ K1 _, v S A6 t8 a7 a
y(N)=z(N);
+ e/ @# q% r$ S
for i=N-1:-1:1
7 N9 X0 h* b' [% X5 S
y(i)=z(i)-u(i)*y(i+1);& x/ I8 H: F* F% m* M4 R, e
end: ^5 G8 h9 M6 z) W; ~1 m2 a
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
& V1 s. ]2 |9 c, r2 y0 |/ H
Y=[af,y,bt]; I v+ g0 z1 Z( d8 K3 Q6 R& E- \
for i=1:N+2
9 H3 f& ^1 J+ [+ {5 X5 y+ B. L3 \
x=a0+(i-1)*h;; _8 S5 x2 y% ?$ Q
zj(i)=1.1392070132*x-0.03920701320/x^2-3*sin(log10(x)/log10(exp(1)))/10-cos(log10(x)/log10(exp(1)))/10;; B3 z1 c+ T7 {
end) L7 ~$ Y$ B9 b
disp('下面两列分别是数值解和近似解');6 Z$ v7 V1 l- X3 m( E, m+ C
re=[Y' zj']