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

[复制链接]

字体大小: 正常放大

2744557306

1189 主题	4 听众	2934 积分

该用户从未签到

电梯直达

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');
( C* _2 R! k5 p) x1 e. d6 a
q=inline('2/x^2');
! X9 n+ W. a0 j4 A/ w
r=inline('sin(log10(x)/log10(exp(1)))/x^2');
( q! ]8 a+ l; W% @, W, R3 p* l
N=9;
1 a/ Q9 {+ ]( Y8 ^2 x
a0=1;b0=2;
& O3 n5 n. ]: n: ~: K
af=1;bt=2;
9 X5 x* s* n8 ^
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
! s( v' x( ]( v$ V$ s
h=(b0-a0)/(N+1);
) @\" G; _) \2 x
x=a0+h;\" G/ ?3 B2 Q3 U2 w. v$ q! y
a(1)=2+h*h*q(x);' ^4 {7 N; W$ U% j4 R5 ?% V4 C
b(1)=-1+(h/2)*p(x);5 l, z& @( |2 ]$ |0 @3 _9 [
d(1)=-h*h*r(x)+(1+(h/2)*p(x))*af;% I; x. e/ g& @9 c* [\" a
for i=2:N-1
$ U. f+ A! `# J/ w# W7 {9 _' y
x=a0+i*h;- G0 J! ^& b\" @% o ?& f
a(i)=2+h*h*q(x);
* E9 ?9 x5 U* j1 M$ f8 q% W
b(i)=-1+(h/2)*p(x);# ] w# x* q3 [% M% N
c(i)=-1-(h/2)*p(x);
: N5 L\" ^) e7 w5 w
d(i)=-h*h*r(x);
v0 Q0 F/ ` m6 ]; l# N& m% @
end
2 P7 W/ v1 m4 u7 u* F
x=b0-h;
+ I9 @4 i2 I4 \+ q5 ?2 T) B
a(N)=2+h*h*q(x);3 C, [+ V& |- u1 g
c(N)=-1-(h/2)*p(x);2 k' G2 F1 p i v\" {8 g; I
d(N)=-h*h*r(x)+(1-(h/2)*p(x))*bt;
- d$ p4 A! @( \( g, c# x
%%%%%%%%%追赶法%%%%%%%%%%%%%%%%%%* o* u! H9 D5 \+ N0 }- r
%y=trisys(c,a,b,d)
5 F# a+ f7 A5 ^/ ^+ b( H+ H
L(1)=a(1);
7 R* M. ^\" p/ N! v5 K
u(1)=b(1)/a(1);
3 J! v( q\" ?0 y9 t% Y\" ?
for i=2:N-1
; O5 m% r) u4 Z1 T1 f$ g. _! Z
L(i)=a(i)-c(i)*u(i-1);0 D$ Q& w) S: d% l
u(i)=b(i)/L(i);7 v8 P% H! o3 v' o: \5 o A
end
- m& ~9 R, @2 f9 Q6 A$ c
L(N)=a(N)-c(N)*u(N-1);
3 o$ l5 J2 D n8 U/ _
z(1)=d(1)/L(1);\" m9 \. L6 P/ ^$ f- L9 F2 k
for i=2:N) W! \, g' J* b8 P
z(i)=(d(i)-c(i)*z(i-1))/L(i);! w; ~. I% B$ @* j, c1 n! `- n
end2 F: Y* q\" m$ w% V+ l/ I
y(N)=z(N);) x# F% W: G+ `; V1 O- z+ W8 ]
for i=N-1:-1:16 P( {' h3 d n# n+ t
y(i)=z(i)-u(i)*y(i+1);
* I1 i' W$ N4 I7 q/ {
end
) L) t- o$ R' V/ J\" u
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%$ b2 a5 Y o# r9 r; a
Y=[af,y,bt];
$ p. ^0 X4 ^\" [% h) U1 m
for i=1:N+27 D, f0 ?, Z% b7 P! O! o
x=a0+(i-1)*h;0 } ~( ]: \. f1 Y1 h; N( v
zj(i)=1.1392070132*x-0.03920701320/x^2-3*sin(log10(x)/log10(exp(1)))/10-cos(log10(x)/log10(exp(1)))/10;
5 \+ r4 M, _# D
end
. W! j' q4 }2 X# Y; l5 o6 `
disp('下面两列分别是数值解和近似解'); N* G2 Z: w, ?6 N! c/ d
re=[Y' zj']