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

[复制链接]

字体大小: 正常放大

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');
8 W$ L' V z8 N% S/ v# D9 T
q=inline('2/x^2');& } ]; U) z/ D3 p\" R' V2 j
r=inline('sin(log10(x)/log10(exp(1)))/x^2');3 \. c, E0 q; e% f2 K( M
N=9;, B/ v4 z n# L* k- i# |2 z
a0=1;b0=2;+ ?0 Q\" n# R. A7 U
af=1;bt=2; R& A1 h% z9 p) J' L: z
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%9 I* ?# E q o/ d5 s: N& E3 F
h=(b0-a0)/(N+1);7 M$ ]9 j7 P% l( l% Q. O
x=a0+h;
q7 y+ c2 t/ X8 v
a(1)=2+h*h*q(x);
' z7 A3 L+ b0 r) a( n
b(1)=-1+(h/2)*p(x);
( v9 s% `' l! x
d(1)=-h*h*r(x)+(1+(h/2)*p(x))*af;5 r2 V! [ F. {\" i# v/ S
for i=2:N-11 W# Z+ F0 X2 M' ?( H$ t, M. D
x=a0+i*h;- t5 P$ E7 W# I
a(i)=2+h*h*q(x);( V; h3 A6 N; Y\" ?; F5 v+ ^8 s
b(i)=-1+(h/2)*p(x);0 k5 o% Y3 H7 a( J F
c(i)=-1-(h/2)*p(x);\" T* g( h' Q. [' t/ y
d(i)=-h*h*r(x);: \- i2 S7 ~5 V8 H
end5 |& V1 R7 b- z/ Z' O
x=b0-h;8 h' B$ c! R0 E& o( X
a(N)=2+h*h*q(x);# F+ H0 W+ _6 x* d2 ~
c(N)=-1-(h/2)*p(x);! O/ v7 n: X& C: |: ]7 F* o
d(N)=-h*h*r(x)+(1-(h/2)*p(x))*bt;
%%%%%%%%%追赶法%%%%%%%%%%%%%%%%%%' C$ J4 L1 r- S' n9 A4 Z$ E
%y=trisys(c,a,b,d)) w' L7 J0 L. B9 W\" s$ B
L(1)=a(1);
2 b6 p2 v) q- g7 N, J
u(1)=b(1)/a(1);. m0 d( }) i4 j+ m7 M) W- k. Y) V
for i=2:N-1
2 g# a- G, U t+ {. h C
L(i)=a(i)-c(i)*u(i-1);7 i! e( H: O8 f1 s\" {3 X
u(i)=b(i)/L(i);
5 G' X! e. [# H5 @
end, U; K; f$ v; b( n2 {; c) A\" |2 j
L(N)=a(N)-c(N)*u(N-1);
3 M( D$ e2 d1 X5 [, z& B
z(1)=d(1)/L(1);- c' g3 g; S! x( }0 r
for i=2:N
1 k4 S\" e& l5 }' e* J0 V' Z\" ^
z(i)=(d(i)-c(i)*z(i-1))/L(i);2 |8 z! b6 |7 p9 Y& G
end
: l; \2 W1 ]* b; I9 |# M& j
y(N)=z(N);
# b% I1 u% `. |
for i=N-1:-1:1
! d' j6 u( c2 _7 x# |! ^
y(i)=z(i)-u(i)*y(i+1);
7 l$ _/ @1 Y5 g- q D. {
end* D: n6 v7 i8 j& w
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
9 Y7 P' {( a. Z* w
Y=[af,y,bt];- z3 E O6 c! Q- y\" S$ `
for i=1:N+2
' v& _- T% _* p% F% l& N) x
x=a0+(i-1)*h;' m' o y6 m' E+ |
zj(i)=1.1392070132*x-0.03920701320/x^2-3*sin(log10(x)/log10(exp(1)))/10-cos(log10(x)/log10(exp(1)))/10;
f7 b7 U; ^& Y
end
6 f9 n) C- a! X1 l5 d
disp('下面两列分别是数值解和近似解');
. m; P! _5 |# x; e- N7 c& I- c& N0 J
re=[Y' zj']