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

[复制链接]

字体大小: 正常放大

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');3 C$ ?4 W! E. o
q=inline('2/x^2');3 Z& S) H a# K! A7 a3 B A
r=inline('sin(log10(x)/log10(exp(1)))/x^2');# ^. U) L( p K$ @# c
N=9;. `; h) _+ v4 z. K) p5 P6 K3 B* {( D5 @
a0=1;b0=2;
. `) N* t$ B! u3 V0 D; c
af=1;bt=2;
+ S' M; C- F( I, D$ L1 }& z7 L
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%! Q2 |) z- D0 u1 g
h=(b0-a0)/(N+1);3 D/ t1 W3 {& j) o5 b! r8 D; j1 i
x=a0+h;( A# J6 X9 k+ H4 Y7 K5 i
a(1)=2+h*h*q(x);0 ~3 e6 t+ B$ r' O, w2 _
b(1)=-1+(h/2)*p(x);; P6 r' Y3 |( a) v1 k+ M$ t
d(1)=-h*h*r(x)+(1+(h/2)*p(x))*af;
/ t% I$ D. K& N
for i=2:N-1
9 Z& U4 E; ?, H7 Y\" `# j
x=a0+i*h;3 Z7 Z4 e3 c( T
a(i)=2+h*h*q(x);4 A( u4 C3 \6 S% R
b(i)=-1+(h/2)*p(x);# p& u% m0 S6 z7 a+ \) s. l
c(i)=-1-(h/2)*p(x);
. W\" o: l6 Z$ W: [+ n! [
d(i)=-h*h*r(x);
9 M! B/ ?( x# ?: d8 V. N' v
end
. h- ^' t% y n+ c- x% Z* G+ F
x=b0-h;' z4 M2 V4 a7 H3 s) U
a(N)=2+h*h*q(x);2 F& f8 Z- N; \- }- a
c(N)=-1-(h/2)*p(x);
) ?0 L3 w\" R+ W c$ D
d(N)=-h*h*r(x)+(1-(h/2)*p(x))*bt;& E4 ?) L\" }# e) \7 e/ b
%%%%%%%%%追赶法%%%%%%%%%%%%%%%%%%% m6 a# j) K( N: j1 a
%y=trisys(c,a,b,d)
4 e. r d8 K7 T( \, m, c
L(1)=a(1);
. V: ~8 j% Z5 o! B' ^7 |
u(1)=b(1)/a(1);\" x3 A) o/ t) d
for i=2:N-17 B0 [( R0 A5 g* {
L(i)=a(i)-c(i)*u(i-1);
$ \% N0 ]9 ?1 j1 l
u(i)=b(i)/L(i);
+ }, `' s3 R6 z* t q G
end
+ d/ b* T& }7 g7 `2 _
L(N)=a(N)-c(N)*u(N-1);
6 b8 V( B* F, @
z(1)=d(1)/L(1);
) Y/ b) y/ i- `* c\" i3 {# k0 _. P
for i=2:N/ t! m( }/ H9 v; p7 h# Z
z(i)=(d(i)-c(i)*z(i-1))/L(i);
$ {2 y( C, u7 ?\" \
end
: V4 ]$ U+ t9 G$ H' c
y(N)=z(N);
- e! f, W% W$ d( p4 v w
for i=N-1:-1:1. h5 S& m& |% K6 s
y(i)=z(i)-u(i)*y(i+1);8 ] b+ ~; g: X# o, L- ]# c( F+ b+ t
end [1 e' b% A X/ Z: z+ Q
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5 u { d5 v# f2 O8 N3 [
Y=[af,y,bt];
, w\" ?5 X0 e2 s( ?8 ~/ c; q+ e
for i=1:N+2! Y& Q$ k* n3 I5 B( i$ ?& S\" n
x=a0+(i-1)*h;
. i! z# x% P\" q$ L' s' |( S( I2 D; g
zj(i)=1.1392070132*x-0.03920701320/x^2-3*sin(log10(x)/log10(exp(1)))/10-cos(log10(x)/log10(exp(1)))/10;2 m% h {1 ]) y9 C2 h1 W
end
& v+ Q1 h- `1 A
disp('下面两列分别是数值解和近似解');
1 b6 G8 D7 P. L( R\" W5 ^* z
re=[Y' zj']