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

[复制链接]

字体大小: 正常放大

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');
q=inline('2/x^2');, S1 S& z2 h4 D
r=inline('sin(log10(x)/log10(exp(1)))/x^2'); s4 E* v' a6 o$ y
N=9;
- h+ z# c0 Y* G' A/ C# v' X% X, ]
a0=1;b0=2;
; c- Q2 W$ v; L% a: V\" }: H1 {
af=1;bt=2;
4 p( g# Y! @! n$ m\" y7 P% ?! {/ H
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%! ~& a) s% l5 B# W
h=(b0-a0)/(N+1);
( s/ o* t( {1 n V
x=a0+h;
: r- A$ C! |$ Q5 o) f5 a7 N
a(1)=2+h*h*q(x);: @. \! ^4 ], d* U7 R4 m* F0 _
b(1)=-1+(h/2)*p(x);
! v9 `8 p( U0 ]: H
d(1)=-h*h*r(x)+(1+(h/2)*p(x))*af;+ |! ^, j( J, O4 G1 v
for i=2:N-18 t' Y# J* p* }+ z
x=a0+i*h;
6 ~3 d6 q( f+ a: W' y2 _, t
a(i)=2+h*h*q(x);' S( l8 @+ {( N+ e
b(i)=-1+(h/2)*p(x);% y% x\" g1 z. m
c(i)=-1-(h/2)*p(x);; N8 k1 E5 ]; ^( ^. P) F: i+ J
d(i)=-h*h*r(x);, b! q% z* M\" m; Q: i6 h8 v$ j
end$ @; W, o6 _1 ^) p. U( T$ }$ H/ ~
x=b0-h;
, ?3 u( n u D/ |; \\" s3 F
a(N)=2+h*h*q(x);
) N* y9 f\" F6 `) D! o3 u& j# A
c(N)=-1-(h/2)*p(x);
9 L7 w1 y* y; ^2 g\" B5 V/ v1 G( M' ^
d(N)=-h*h*r(x)+(1-(h/2)*p(x))*bt;
/ R5 X\" t3 S# t1 N3 ]- W4 o
%%%%%%%%%追赶法%%%%%%%%%%%%%%%%%%
- s2 K1 l3 `' c7 z; W5 q4 c3 U
%y=trisys(c,a,b,d)
8 S7 j9 M9 M# d4 Y; c
L(1)=a(1);2 _$ w8 F4 n# a& t9 n
u(1)=b(1)/a(1);* ] T+ [: g/ l( I( ~0 b2 F. K
for i=2:N-1; P9 {) B$ V6 V$ j0 e4 E% O2 [2 b
L(i)=a(i)-c(i)*u(i-1);1 \\" f, K/ H1 y6 O+ }0 \
u(i)=b(i)/L(i);
. t2 b1 i( V! }6 I
end+ B) {* j6 ?: S6 `/ J
L(N)=a(N)-c(N)*u(N-1);
8 e6 K; ~3 G0 E
z(1)=d(1)/L(1);
- W5 d) j: G; j( M( h
for i=2:N1 R1 u2 p: ]* |; a$ K- H% D
z(i)=(d(i)-c(i)*z(i-1))/L(i);$ f3 C\" ^/ ~7 ~, r' `# t7 {
end% a$ V) C; \! T) L. ~; w# g( q$ I
y(N)=z(N);) P- ~2 j4 U. L! S$ `9 H! G
for i=N-1:-1:1$ w# D! @( t9 T6 U0 |
y(i)=z(i)-u(i)*y(i+1);
Q( p% n8 [, |' L
end' r& J+ U0 h- x\" n! h% ]- A) O
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5 ?; U8 d$ ]6 K! ^: u9 c! g/ F
Y=[af,y,bt];
4 J$ x! R\" f- B% S3 b: T& }9 m
for i=1:N+2
. j' ~) p7 v! i* S4 r
x=a0+(i-1)*h;
( k! j% d* w H2 }2 \3 b3 \9 m
zj(i)=1.1392070132*x-0.03920701320/x^2-3*sin(log10(x)/log10(exp(1)))/10-cos(log10(x)/log10(exp(1)))/10;
- D! ^$ C c% N/ K! _0 c: @& N6 E3 v
end& R' ^; r& E \+ l+ @* U
disp('下面两列分别是数值解和近似解');
5 r7 N( ]9 ?& q7 Q
re=[Y' zj']