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

[复制链接]

字体大小: 正常放大

2744557306

1188 主题	4 听众	2931 积分

该用户从未签到

电梯直达

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');1 ]0 x; d1 b\" W5 ?5 q+ x2 w* I
q=inline('2/x^2');
0 H5 J\" U+ L$ C: n, p
r=inline('sin(log10(x)/log10(exp(1)))/x^2');& w9 y% h* S- x! u, A& n+ t
N=9;/ n8 c$ E- ]( r3 u6 |& _& k2 l
a0=1;b0=2;
/ @0 Q( p) K% w
af=1;bt=2;
7 I) p i, t- p
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%; D( R9 k, [4 f. k7 _; m8 P% B
h=(b0-a0)/(N+1);
& z5 [3 R2 G1 B1 r1 \
x=a0+h;: x. Y9 _! K- a9 Z
a(1)=2+h*h*q(x);
) ~- h4 _# h' o# u$ ^
b(1)=-1+(h/2)*p(x);
, b g# G: ]% m$ d, H
d(1)=-h*h*r(x)+(1+(h/2)*p(x))*af;
; X/ T4 r$ S4 P9 v! o7 @7 O# P. j
for i=2:N-1
- e3 b2 ]* T- t/ V& f8 s
x=a0+i*h;( j% y/ Q/ y1 o\" A7 i
a(i)=2+h*h*q(x);
+ ?+ z! y* |% c$ J, U
b(i)=-1+(h/2)*p(x);
. E& Y% w7 j* \6 [, w
c(i)=-1-(h/2)*p(x);
P6 h/ P8 Y& ?5 \9 K' N6 W
d(i)=-h*h*r(x);
; \4 y @4 ^& b1 K# J+ w9 d9 w# h
end
/ h. `3 r# A8 u
x=b0-h;3 _! D/ r\" u2 A
a(N)=2+h*h*q(x);
c(N)=-1-(h/2)*p(x);3 {, w9 @9 R! J4 i
d(N)=-h*h*r(x)+(1-(h/2)*p(x))*bt;/ q4 c, F$ Q& y5 V/ @* O2 e
%%%%%%%%%追赶法%%%%%%%%%%%%%%%%%%
; {9 s; _( j# j0 p. `
%y=trisys(c,a,b,d)
5 h' s+ q- W* h! H
L(1)=a(1);
' u) G' @% I$ X- s
u(1)=b(1)/a(1);
* T- Q* b. N0 a7 _% `* j/ F @8 ?
for i=2:N-1+ |: J$ _1 u! W, M. N% T% q: s0 Y `& F
L(i)=a(i)-c(i)*u(i-1);
0 \: i/ x: i& ~ e
u(i)=b(i)/L(i);
& {5 y' j) c: d
end
# d1 G. l5 b# r- U1 C
L(N)=a(N)-c(N)*u(N-1);
7 v5 j: g d7 o
z(1)=d(1)/L(1);
9 R% u1 |1 X* c. I' n
for i=2:N
/ n& L4 k' i& b( i
z(i)=(d(i)-c(i)*z(i-1))/L(i);% |* n% y% N# d& Q
end
, I8 c1 \. x) i# }' s
y(N)=z(N);1 R3 ^- ]8 L; ^) b1 W5 X* l$ |9 u
for i=N-1:-1:1
; c! G+ j7 q9 B1 J6 ~1 w8 W
y(i)=z(i)-u(i)*y(i+1);
! s/ c0 T6 v% l: a6 N8 z9 T2 w) T, [
end, ?8 r6 K2 a, z3 G$ H
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ b1 t* g M7 P3 i! ]
Y=[af,y,bt];
; C# J3 \9 `* d+ R. s( t* u& O( R
for i=1:N+27 @0 A$ f9 u3 G1 P7 C1 K* N
x=a0+(i-1)*h;+ h) m# X' y- z
zj(i)=1.1392070132*x-0.03920701320/x^2-3*sin(log10(x)/log10(exp(1)))/10-cos(log10(x)/log10(exp(1)))/10; x1 ~. W% y6 S! z+ s: p! F. }+ e7 S\" C' h
end
4 U\" O: n# d1 j( e9 J6 f
disp('下面两列分别是数值解和近似解');( e; s( q9 x- b- ]- q( L
re=[Y' zj']