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

[复制链接]

字体大小: 正常放大

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');! z( S\" d3 t: @% E7 \& d I# `
q=inline('2/x^2');
- {) d; O$ I. b4 S8 e
r=inline('sin(log10(x)/log10(exp(1)))/x^2');
6 ], I! _% h\" Y3 B/ d+ \6 Y2 L
N=9;2 d3 y) F; Y2 Q8 [! x
a0=1;b0=2;6 U% p% S0 N; t$ c, z3 g, v
af=1;bt=2;
1 V6 i( c( T' z2 b S( R' l
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%0 b2 N: D. K4 \) k* t, }$ C9 W
h=(b0-a0)/(N+1);\" I8 S3 J e+ k u
x=a0+h;
) ^5 g. k7 y% ?: ]+ y0 q
a(1)=2+h*h*q(x);
3 K\" N: ~5 Q$ k+ _
b(1)=-1+(h/2)*p(x);
8 B, \$ {* \2 X' [
d(1)=-h*h*r(x)+(1+(h/2)*p(x))*af;
) ~1 M( l( l& p
for i=2:N-1
- z- w. g) H: c
x=a0+i*h;
# e& F O# D. l
a(i)=2+h*h*q(x);
0 Z0 W2 C) c2 }: _- M) B
b(i)=-1+(h/2)*p(x);
: P9 A$ [! q( A y2 S
c(i)=-1-(h/2)*p(x);: j$ n0 t5 B% I! ]8 p
d(i)=-h*h*r(x);7 Z# M\" r9 O. V6 I+ F
end
2 q7 [. ?: S0 U2 W7 t K8 z
x=b0-h;! P! l8 n- F7 ^) I) L
a(N)=2+h*h*q(x);3 M& h9 J! P7 l
c(N)=-1-(h/2)*p(x);
/ S0 K- Y, d, w8 k6 G: A- P
d(N)=-h*h*r(x)+(1-(h/2)*p(x))*bt;- U2 v' f! M7 B
%%%%%%%%%追赶法%%%%%%%%%%%%%%%%%% W: d# e! I9 A2 o
%y=trisys(c,a,b,d)
7 e0 ~5 p0 S\" V
L(1)=a(1);9 t0 j\" [7 i1 R1 O1 c# X
u(1)=b(1)/a(1); f( R. _4 Q0 ?% K7 @
for i=2:N-1
1 C/ I- N; Q* B5 J
L(i)=a(i)-c(i)*u(i-1);, H* U4 w2 }# s- o4 Y8 `3 J4 L
u(i)=b(i)/L(i);
& r3 B/ k& Y3 S0 u) T
end+ T9 \' l\" B8 ]1 v! A8 H0 v, g
L(N)=a(N)-c(N)*u(N-1);\" \; c) q: k. R2 W- k Z9 S4 Z; N9 m9 _
z(1)=d(1)/L(1);
7 @+ O$ `\" e& D: \ O) A
for i=2:N1 R3 o; C2 t0 F3 @( f9 ]6 @
z(i)=(d(i)-c(i)*z(i-1))/L(i);
. d, y\" J4 c+ T, q
end$ r\" U\" U; B U9 r
y(N)=z(N);& W2 Z7 s- ^1 u# Y/ R! ^
for i=N-1:-1:1; T1 I9 W' G' |) q, G
y(i)=z(i)-u(i)*y(i+1);
- v, R1 y2 z- P! P, V: g
end: S* K% S$ s' I' z0 B+ O
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
( ]$ V% z* D! w% [! Z) m8 c
Y=[af,y,bt];. b* w7 ^- N6 t& S' W. H) }- W+ z* D
for i=1:N+2
! Z* E& x' x4 O! s' }% @0 T5 t
x=a0+(i-1)*h;\" I! s\" ?3 b: l4 Q/ t7 ~2 f/ R) ]* t9 d
zj(i)=1.1392070132*x-0.03920701320/x^2-3*sin(log10(x)/log10(exp(1)))/10-cos(log10(x)/log10(exp(1)))/10;- r) W3 R' p T/ r& `
end
$ p4 k- Q7 u* T, z* L0 y8 C$ r
disp('下面两列分别是数值解和近似解');
) C- q- v. D( [& i5 g0 n; C
re=[Y' zj']