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

[复制链接]

字体大小: 正常放大

2744557306

1186 主题	4 听众	2922 积分

该用户从未签到

电梯直达

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');+ f$ t/ |% }3 O- J! d# p\" t4 h
q=inline('2/x^2');' \- S u: ?( x: R
r=inline('sin(log10(x)/log10(exp(1)))/x^2');' V6 i {1 P0 ^5 s9 K4 G( Q3 p
N=9;
. k5 Y( `7 V8 V2 X' b
a0=1;b0=2;
4 K) g' ]6 N2 R* v( k
af=1;bt=2;
3 _8 y0 Q* w' g' W4 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%; {9 R R- K h2 b: q# N
h=(b0-a0)/(N+1);
$ v! b1 T, L) @6 J# H. @5 t
x=a0+h;4 t1 o; D0 A* j& u2 U8 [
a(1)=2+h*h*q(x);3 @* U6 |0 N) a) D( h' V
b(1)=-1+(h/2)*p(x);: N\" H. N) @; D
d(1)=-h*h*r(x)+(1+(h/2)*p(x))*af;
& t R+ M, C$ n9 |# [! s
for i=2:N-1
; S; t- V0 y* o0 E. X7 [2 W
x=a0+i*h;
& l, {1 K, w: O$ c% z9 n
a(i)=2+h*h*q(x);; T: a- m- {, j2 `
b(i)=-1+(h/2)*p(x);6 r1 \5 g1 E$ a7 m' f
c(i)=-1-(h/2)*p(x);
! p5 V$ S2 W! }
d(i)=-h*h*r(x);
L# r, O. }! x
end
% m$ Z2 F7 m( i( A% L8 _) g# v
x=b0-h;0 G8 ]. W' d m/ D$ U9 \
a(N)=2+h*h*q(x);6 u: a' G; h4 z9 a2 Z
c(N)=-1-(h/2)*p(x);; w9 u# R% n3 b* R* f: A
d(N)=-h*h*r(x)+(1-(h/2)*p(x))*bt;
6 W. j7 W c# W8 O1 R! |2 ?
%%%%%%%%%追赶法%%%%%%%%%%%%%%%%%%; j/ x4 z* @6 C+ `4 ~ J) u
%y=trisys(c,a,b,d)9 Y! P$ t0 ], C
L(1)=a(1);$ [1 d! {\" E$ i9 D: p
u(1)=b(1)/a(1);
: f* R/ d ~3 Y' p/ ]; o0 ` E# }
for i=2:N-1# @- V; G5 \/ `6 F
L(i)=a(i)-c(i)*u(i-1);1 I1 }+ P# {9 e8 H
u(i)=b(i)/L(i);
; C\" E4 u S, ?3 l
end- X6 h& \+ P+ U; q6 T
L(N)=a(N)-c(N)*u(N-1);
& C6 M8 w0 e0 d! T- ?9 z6 Z
z(1)=d(1)/L(1);' {1 A0 b$ g# H; E, T4 U
for i=2:N6 t) m2 i. K3 v: ?; _; ^
z(i)=(d(i)-c(i)*z(i-1))/L(i);) _6 Y9 D1 U$ K
end
$ ^' ~4 Y Q9 G
y(N)=z(N);
5 v: T0 H1 ~/ Z* s3 x
for i=N-1:-1:1
- I+ b8 Y4 z# Q+ R
y(i)=z(i)-u(i)*y(i+1);+ F+ k: y. H& c+ ]
end, h! B' v! j0 i2 X( g( B& h
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
! p% _: {( p O3 [1 F
Y=[af,y,bt];0 Y% H5 k1 C7 `3 x\" @
for i=1:N+2
* C: a. |9 ?9 g' d; T
x=a0+(i-1)*h; h6 g q: i, G/ w$ \3 u
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$ j+ o/ |6 g
end6 y2 o3 @' u4 o J
disp('下面两列分别是数值解和近似解');
+ J; B1 d5 t4 L; }/ m- k
re=[Y' zj']