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

[复制链接]

字体大小: 正常放大

2744557306

1192 主题	4 听众	2946 积分

该用户从未签到

电梯直达

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');0 ~. K9 ]5 H% `\" L5 d
q=inline('2/x^2');
, b, \, Y\" q& W# r
r=inline('sin(log10(x)/log10(exp(1)))/x^2');3 _6 i: C/ ]; V+ d; u. s
N=9;1 B8 Z\" T% i; n/ m1 z
a0=1;b0=2;
# l( b0 R$ g: F9 @! {3 G
af=1;bt=2;
) `8 e2 M/ v' I
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# @! R( h H0 U/ I
h=(b0-a0)/(N+1);
+ L% \/ ]: B' R7 g6 J4 o( f7 o
x=a0+h;. e# L4 _; [7 i\" O
a(1)=2+h*h*q(x);9 a& {1 h& Y7 E: h
b(1)=-1+(h/2)*p(x);
! M$ [8 C1 Z9 c% A
d(1)=-h*h*r(x)+(1+(h/2)*p(x))*af;
& x O6 M/ u! @3 ?
for i=2:N-10 ~# h$ W0 {0 _' n7 ]4 I7 @2 w
x=a0+i*h;
; F, Q$ g- L+ b6 @
a(i)=2+h*h*q(x);$ \1 r1 P/ ^6 \) Z: }* l( A9 s& h
b(i)=-1+(h/2)*p(x);
* ~% z* O C1 s4 ^
c(i)=-1-(h/2)*p(x);. O0 |& P; {( _5 B% [9 x& w
d(i)=-h*h*r(x);
4 D1 y8 V: ~8 j
end, o- Q4 o8 _0 |8 c9 N# S
x=b0-h;& _& t, }1 ?; y& C j A
a(N)=2+h*h*q(x);0 C/ o3 t2 S/ \0 q% `
c(N)=-1-(h/2)*p(x);
9 W# T7 b3 [! ^# A+ a5 W8 D2 d$ d1 d5 R
d(N)=-h*h*r(x)+(1-(h/2)*p(x))*bt;3 _$ X0 x. z4 {0 d& M3 L$ S+ y
%%%%%%%%%追赶法%%%%%%%%%%%%%%%%%%
! I+ X! [\" T! @/ D
%y=trisys(c,a,b,d)) B0 `( ]! _3 E
L(1)=a(1);
. @. @/ B: M2 }. D6 l1 b
u(1)=b(1)/a(1);# b5 s& m0 I. f: a
for i=2:N-16 u P/ X7 q; B
L(i)=a(i)-c(i)*u(i-1);
0 F) `4 p8 W$ V. y5 Z* E$ Q
u(i)=b(i)/L(i);
! r4 Y5 @' C: ]3 J
end9 d. x$ w; r+ w ~, z
L(N)=a(N)-c(N)*u(N-1);
7 k, ?; C' K( t4 Z; D) x
z(1)=d(1)/L(1);
% z: R* C6 a' k
for i=2:N
! F7 P$ ^2 I0 j7 O: T# m
z(i)=(d(i)-c(i)*z(i-1))/L(i);
& ]5 J5 a* O) Z. }$ f, e
end
4 s& j' {) ]\" e i. o\" `$ m1 L3 f
y(N)=z(N);
; G% T6 d- ~: j. \& q: Q\" ]
for i=N-1:-1:1
. i5 N2 q9 V4 K1 K5 \ U6 o; K
y(i)=z(i)-u(i)*y(i+1);, J* s) ]- f+ Z- C+ D
end9 a+ M# U4 L: P3 g4 G/ Y# p8 \
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%2 N8 ^2 O, a$ b A- ?6 p
Y=[af,y,bt];
# A' ~ `2 g6 A, s. D% [0 z# F
for i=1:N+2
/ G+ \, x1 b3 N7 G
x=a0+(i-1)*h;* p# H! Q7 \( f2 }! u2 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;! h! D8 ^% G# o2 S# r- T. C8 l- b2 r
end
! i- `$ N3 @8 m, o& E# Y6 Y
disp('下面两列分别是数值解和近似解');; n& ?7 r; S0 M7 |# J3 L! [' v' u
re=[Y' zj']