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

[复制链接]

字体大小: 正常放大

2744557306

1184 主题	4 听众	2916 积分

该用户从未签到

电梯直达

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');
, d) P4 R! E# h/ W' \0 M
q=inline('2/x^2');& ?1 H- P& m$ e( [% ~5 { J, Y
r=inline('sin(log10(x)/log10(exp(1)))/x^2');8 I) y( a; N. K/ w\" J( n o
N=9;5 q; h5 N7 d& `4 @$ p% N
a0=1;b0=2;4 z4 d4 n @ @% @
af=1;bt=2;
1 E( D: m8 Q8 Y* a6 j\" Y
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%# J5 S: m+ F7 c1 Z
h=(b0-a0)/(N+1);
4 C. v# X5 R, _' z
x=a0+h;
) b4 N) Q; {- | `2 e: {
a(1)=2+h*h*q(x);( a# @! P+ s$ h$ z( K1 o1 D
b(1)=-1+(h/2)*p(x);
) W+ V( H' C0 Z8 f7 C& y. i
d(1)=-h*h*r(x)+(1+(h/2)*p(x))*af;
\" D' G+ D- q2 c
for i=2:N-1
# }9 @. M( B: o% E6 o! F
x=a0+i*h;8 Z$ W9 Y5 r9 M\" H l7 J; k6 E
a(i)=2+h*h*q(x);
# P$ w7 E. Q1 D' y* {' K- x
b(i)=-1+(h/2)*p(x); ?% V& N0 O' K3 N# A& c; I$ p4 `) p1 j+ B
c(i)=-1-(h/2)*p(x);
& w0 d( Q$ |8 a, Y4 a' |* o& Y
d(i)=-h*h*r(x);$ Y% U) }; x: l- {
end2 L\" B4 a A2 q5 E6 P
x=b0-h;
$ g/ T9 f3 S\" l- b* S
a(N)=2+h*h*q(x); @& X7 N9 U; r. o\" T' w: Q) s
c(N)=-1-(h/2)*p(x);\" r) m' h& R\" @' I0 o7 x
d(N)=-h*h*r(x)+(1-(h/2)*p(x))*bt;
+ l) Z9 ^( p0 Y: r+ S\" R3 \
%%%%%%%%%追赶法%%%%%%%%%%%%%%%%%%
5 d4 P3 l! c4 x- V
%y=trisys(c,a,b,d)' P& f& ]! }, }2 N: K5 [
L(1)=a(1);+ P5 c) p7 v2 r7 a* C% W
u(1)=b(1)/a(1);
% r* Q( |/ r! [: N
for i=2:N-1
5 J0 y4 j, {3 `! c
L(i)=a(i)-c(i)*u(i-1);
h- g. m7 {\" K% `
u(i)=b(i)/L(i);
}$ R- Y' T) x$ y9 r- Q* W, u\" u
end
0 o$ e' K( G- X$ Q
L(N)=a(N)-c(N)*u(N-1);
% l3 r3 h0 N# X I; s- H8 W
z(1)=d(1)/L(1);
7 b; }( g: y/ _2 W& N% d
for i=2:N
\" g9 B! P# v3 F( `* A
z(i)=(d(i)-c(i)*z(i-1))/L(i); ]. _( a$ D# X& j6 N5 l6 p
end9 s\" T9 _4 z# F: S* x
y(N)=z(N);
: s- C; |8 ~1 ]\" x, `% @
for i=N-1:-1:14 ^: p4 v4 D, \! J) d( B
y(i)=z(i)-u(i)*y(i+1);* I( F, l4 [$ v( Y) e U
end
0 M- I o! |! M; N, | f
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%. `+ E! ~' [; y\" B& Q- H
Y=[af,y,bt];3 U1 M0 ^1 B/ A `0 U; p
for i=1:N+2
r/ E6 v/ e: p! j$ j8 K+ K( b
x=a0+(i-1)*h;& l5 K. x3 v0 X, j, {8 k7 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;( c( Q\" {( X( [& g\" X( z. P
end$ r5 G9 h v, o\" j+ g
disp('下面两列分别是数值解和近似解');! F& O* c3 W& \' K6 _2 G( b
re=[Y' zj']