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

[复制链接]

字体大小: 正常放大

2744557306

1198 主题	4 听众	2978 积分

该用户从未签到

电梯直达

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 g6 s8 S( G* W$ \8 n% g
q=inline('2/x^2'); n, j# v* x2 r. D
r=inline('sin(log10(x)/log10(exp(1)))/x^2');) d4 L+ ^, m8 I5 }) V. c% {7 W
N=9;
0 R+ z ^' {$ O, l
a0=1;b0=2;
$ V( H% N% j @0 S) s2 @, \$ R4 v
af=1;bt=2;7 r N/ o& \4 Y\" x, x$ J2 P
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\" \% V( R: H# l( I. [
h=(b0-a0)/(N+1);\" O% Y6 v7 H! `% L# {. k6 N4 }# S\" E
x=a0+h;' h$ N8 ~7 }# G# M0 ^# R
a(1)=2+h*h*q(x);
6 z/ s\" g6 H( A8 x! u
b(1)=-1+(h/2)*p(x);
. r7 Y( W3 p5 [* ]2 U+ e4 B- v
d(1)=-h*h*r(x)+(1+(h/2)*p(x))*af;0 I2 H' O) |- _( ?
for i=2:N-1$ D, t8 r1 T( W3 e7 j2 h; r8 ?
x=a0+i*h;& f3 Q, Y( G# E, @# j8 f
a(i)=2+h*h*q(x);
, ^* f( Q5 X1 d N% O\" [. B
b(i)=-1+(h/2)*p(x);1 o2 i. B& k0 q/ f& K% Q
c(i)=-1-(h/2)*p(x);
( X; G\" Y# r* r; V& p7 s1 Z+ g
d(i)=-h*h*r(x);
4 Q: G: |. t. F2 ~
end; N# M: Z3 W& x+ E/ o# a
x=b0-h;1 x2 a4 n5 X. z4 v$ f2 F. c l
a(N)=2+h*h*q(x);
9 [) r T' m& g/ \* @\" @, \. k1 z
c(N)=-1-(h/2)*p(x);3 Y+ [7 T- K; m! F& o0 \# q
d(N)=-h*h*r(x)+(1-(h/2)*p(x))*bt;
3 `( V8 p. Q) X8 Z
%%%%%%%%%追赶法%%%%%%%%%%%%%%%%%%( b4 l) {% y1 ^, o
%y=trisys(c,a,b,d)- C! h O+ \: @1 y
L(1)=a(1);1 [! |3 x, N- b s6 ^! `9 ^) z
u(1)=b(1)/a(1);
+ I$ R6 O0 c9 I! i
for i=2:N-1
! R/ s& d( X2 g2 A' T: C
L(i)=a(i)-c(i)*u(i-1);
\, L2 ]5 q k) V
u(i)=b(i)/L(i);
/ Q- |! D4 d1 y) A6 O( V
end, f, |5 V5 a0 x0 E
L(N)=a(N)-c(N)*u(N-1);4 w; o F( p\" u+ j7 H\" H
z(1)=d(1)/L(1);& ^6 j- }& y& ]$ \) l
for i=2:N2 H$ F& A7 y. j7 x' H
z(i)=(d(i)-c(i)*z(i-1))/L(i);
4 z* S9 |' u, O, R5 o5 s
end
2 F P* A1 p1 c5 c7 k @; `, Q' Y
y(N)=z(N);7 s5 b: n! X8 d\" i% e
for i=N-1:-1:1
$ ^* s# @! D8 F* h
y(i)=z(i)-u(i)*y(i+1);
0 g, N |+ O+ Z( P
end2 S1 @' N: j+ c9 e: w
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4 ]7 G/ d; L% o- t9 [7 P
Y=[af,y,bt];
- x$ \* x0 C! M6 e. j& p
for i=1:N+22 c+ T- ^# k6 G
x=a0+(i-1)*h;
\" Q4 e6 T6 Q9 s8 r5 _: F. B
zj(i)=1.1392070132*x-0.03920701320/x^2-3*sin(log10(x)/log10(exp(1)))/10-cos(log10(x)/log10(exp(1)))/10;7 X! ]( V& L: v) ]! q0 E
end
5 P# e$ C0 { W: H; `* e\" L
disp('下面两列分别是数值解和近似解');8 \+ D6 |$ h$ L9 @+ S4 G& Z& `3 Q
re=[Y' zj']