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

[复制链接]

字体大小: 正常放大

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');
( ?$ `4 Z& D\" m! U( @: O
q=inline('2/x^2');7 K. v. }9 o8 [2 A- v
r=inline('sin(log10(x)/log10(exp(1)))/x^2');6 U( L/ D9 K9 c( [. w
N=9;\" K% X# F7 b! e7 P
a0=1;b0=2;
; H- J5 T ^* B
af=1;bt=2;- K k; F3 f5 i/ f2 k9 Q
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%1 h- T, d; q: I; L7 i, d# U
h=(b0-a0)/(N+1);
, O: r' q( c6 U
x=a0+h;* a; i; \' P& D- b: f
a(1)=2+h*h*q(x);
\" n I7 ]/ V: n9 Z/ E
b(1)=-1+(h/2)*p(x);
. u+ s7 L, E5 ~0 V7 x, j
d(1)=-h*h*r(x)+(1+(h/2)*p(x))*af;
- M' \5 [: }; j/ {' J4 D. H
for i=2:N-1
! m2 L7 x\" S! D1 h\" F( R
x=a0+i*h;
# j( B\" T( }. ]7 B
a(i)=2+h*h*q(x);
$ ]2 q {+ L& c& e) Q* j% X) m
b(i)=-1+(h/2)*p(x);( I k% d4 x7 g2 n4 M
c(i)=-1-(h/2)*p(x);
' r2 Z: Z! |, v! U
d(i)=-h*h*r(x);
; e3 Y7 y9 k) ^& i5 d' M
end; C! f1 m) G+ e- T
x=b0-h;\" v0 l3 Z! _+ }: a7 m- m; j R
a(N)=2+h*h*q(x);
$ k8 Y1 f% j) z& I5 \0 Q
c(N)=-1-(h/2)*p(x);
5 z+ T# R$ d# |& M& W4 H, S) g\" V( v
d(N)=-h*h*r(x)+(1-(h/2)*p(x))*bt;
% f* Q$ \( W7 Q2 G
%%%%%%%%%追赶法%%%%%%%%%%%%%%%%%%\" n3 y8 y+ U# q Q- A# t1 Y, J
%y=trisys(c,a,b,d)
/ U, \, M$ [1 v\" h* J: l
L(1)=a(1);
z; @- C: C6 }) D
u(1)=b(1)/a(1);- C% @% i4 w4 u6 F. |\" `+ y
for i=2:N-11 E\" z7 G7 L: h6 W
L(i)=a(i)-c(i)*u(i-1);
- p; ^4 s+ k5 c, E% R7 W2 v6 L
u(i)=b(i)/L(i);
9 e+ j: {- k3 h9 }# S2 ~
end
# x _4 n) i& \: T, u
L(N)=a(N)-c(N)*u(N-1);
9 ~3 A: H. ^2 l, j. d8 }
z(1)=d(1)/L(1);% W0 F, n) B% J0 r# w\" G( [
for i=2:N5 O+ F. t% [! X! w6 z: b- a
z(i)=(d(i)-c(i)*z(i-1))/L(i);
% r2 p# a/ \5 I4 |
end7 u) A! h5 Y d$ x6 w* g1 u1 q4 ^
y(N)=z(N);
\" f' w+ j z q7 ?% x
for i=N-1:-1:1 w) ?, }: W# u% ?2 i\" S; W
y(i)=z(i)-u(i)*y(i+1);' U; d4 f1 ]1 x; u3 |/ k% x
end q9 o8 C7 G1 h) U0 t$ ?; v\" v: l
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%! M V6 z9 ]7 y( L: H# i6 v
Y=[af,y,bt]; W% D, ^7 h0 ~# E! v' K1 l\" q+ b3 K
for i=1:N+2% r. t\" v/ R( K# ]& ~
x=a0+(i-1)*h;
0 b0 r6 g4 A. [* t! l
zj(i)=1.1392070132*x-0.03920701320/x^2-3*sin(log10(x)/log10(exp(1)))/10-cos(log10(x)/log10(exp(1)))/10;' u4 ]( A5 t. c, X6 p+ t) U K8 @
end* @\" s% [+ O6 W8 }) l4 s' C
disp('下面两列分别是数值解和近似解');6 ~1 H/ t; a4 _1 {/ @
re=[Y' zj']