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

[复制链接]

字体大小: 正常放大

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');5 g3 Y' r5 Q3 h
q=inline('2/x^2');/ s3 ?, F4 O8 y# h
r=inline('sin(log10(x)/log10(exp(1)))/x^2');
! a5 k }& y3 D8 C4 K
N=9;
1 D( ^# y8 ]( t2 X
a0=1;b0=2;- Z0 C+ M4 i\" B' a
af=1;bt=2;
6 v( D. a. Y; W
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
( t, G J& J0 Q
h=(b0-a0)/(N+1);
, q& |2 U! M' @6 K$ q
x=a0+h;
3 _, ^4 k& E# e3 Y: O( |: t M
a(1)=2+h*h*q(x);8 B' [9 O, P) i- I2 N( e
b(1)=-1+(h/2)*p(x);
* { x& p3 ^* N1 d! z0 }0 ~! b
d(1)=-h*h*r(x)+(1+(h/2)*p(x))*af;
0 S3 T+ M8 O* X. S& p
for i=2:N-14 `/ M) ^. l. h4 |! ]% P4 n
x=a0+i*h;
' x8 u' H3 o. [9 T% m! {
a(i)=2+h*h*q(x);! K* w\" C3 y! Z
b(i)=-1+(h/2)*p(x);
6 | ~9 X: B+ _/ `% |% h. l( X
c(i)=-1-(h/2)*p(x);6 Q0 l1 v+ [& [ `
d(i)=-h*h*r(x);' M- N4 o6 }, z- t- d' m6 Y
end) e# m g' e* i# K
x=b0-h;: W% V4 G% B* G% z. l9 j) I3 R
a(N)=2+h*h*q(x);! Z2 }5 [# y. n# U
c(N)=-1-(h/2)*p(x);% _. ~$ m/ x4 N9 [! \8 O
d(N)=-h*h*r(x)+(1-(h/2)*p(x))*bt;
) n T @5 x ]. m) P* q
%%%%%%%%%追赶法%%%%%%%%%%%%%%%%%%
9 Z( ^- p' N- t, T) @4 |
%y=trisys(c,a,b,d)- W0 N( v2 c6 a) G, A' b% Y; q
L(1)=a(1);
! c& o7 V# h6 [1 T* P
u(1)=b(1)/a(1);; V$ j- S2 a/ u. e5 l. r0 x* W
for i=2:N-1
. C9 O6 ~3 C/ p2 |
L(i)=a(i)-c(i)*u(i-1);1 {0 B2 h8 @( }+ A9 j: ?\" @
u(i)=b(i)/L(i);
\" A- R( p2 p- M6 p8 {7 V: Z# f7 Q
end
% G) G5 E4 A: e0 Q {
L(N)=a(N)-c(N)*u(N-1);
; i4 J* C2 e3 u
z(1)=d(1)/L(1);
# e% @. c& G& _% K& r& b: A% D
for i=2:N( c6 s x\" l, _; T5 r
z(i)=(d(i)-c(i)*z(i-1))/L(i);
, Y2 s/ `4 C. b, B+ W
end
1 D* u8 p* W( e
y(N)=z(N);- `7 r# o! }! Y
for i=N-1:-1:1# w5 y; C& _: U% O# G
y(i)=z(i)-u(i)*y(i+1);
% U+ k( I4 n! N
end/ s# h2 M c t6 V' s. i8 c* ?
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
* R4 M& J$ o2 D. z# R. p: _
Y=[af,y,bt];
, J) W. X0 [& T, l3 V# {) P
for i=1:N+2. [. p. f8 s/ k
x=a0+(i-1)*h;
' B& Q( G9 i: T( V/ t {
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 j) p7 I$ R g w3 I) s) f, r# Z
end+ L+ ^! x' \. X% B
disp('下面两列分别是数值解和近似解');; V' @0 [# N$ g: O
re=[Y' zj']