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

[复制链接]

字体大小: 正常放大

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');
: Z3 y6 K( }' A S: Z\" D0 \
q=inline('2/x^2');
* A; b) N+ F, d: D7 q& [8 W& n- S
r=inline('sin(log10(x)/log10(exp(1)))/x^2');# n2 n: ^% j9 z3 x, c1 P
N=9;1 x2 T/ M\" r; b4 N/ r
a0=1;b0=2;& q3 ~: Q. ]8 b+ G1 K' ^( X+ Y
af=1;bt=2;; [/ x4 q$ r8 J/ W
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2 K6 q k* b$ w4 s
h=(b0-a0)/(N+1);. ?+ k4 L/ g J# W& n2 N/ J2 g
x=a0+h;& {0 g: J\" |9 Y
a(1)=2+h*h*q(x);
) e3 E# E( r5 b
b(1)=-1+(h/2)*p(x);8 k1 g8 l5 m9 }
d(1)=-h*h*r(x)+(1+(h/2)*p(x))*af;
+ L; s D4 r1 R/ [& {
for i=2:N-1+ ?/ M$ I; Z: C2 g$ ~- u
x=a0+i*h;
( h% P) n, }. t* O, @# y
a(i)=2+h*h*q(x);- x9 _1 @+ R' E# X/ i1 x; m7 a
b(i)=-1+(h/2)*p(x);: g: p4 D. e5 O: D/ B
c(i)=-1-(h/2)*p(x);
: ]1 `7 }: \6 D
d(i)=-h*h*r(x);
$ F, M5 O& @% X9 P' P7 Y$ c
end, M3 _* z- {5 y( _0 e3 W% r4 v6 F
x=b0-h;
8 s2 A6 D7 C4 x! T3 [
a(N)=2+h*h*q(x);& A! I\" k9 i5 k3 x* |% m
c(N)=-1-(h/2)*p(x);9 O* W* @; M5 R- r$ M1 X
d(N)=-h*h*r(x)+(1-(h/2)*p(x))*bt;3 A. t) f2 M- ?9 V. q' \2 C
%%%%%%%%%追赶法%%%%%%%%%%%%%%%%%%; J4 r3 G8 f. Y( n* W
%y=trisys(c,a,b,d)
0 |+ | l\" M1 ]( j
L(1)=a(1);
( K1 D. D6 N5 I# A; C# o1 K
u(1)=b(1)/a(1);* _1 H% t' k$ L- q3 K
for i=2:N-1
; @3 Z7 X6 N8 x2 @\" I
L(i)=a(i)-c(i)*u(i-1);
, @0 X/ r9 v& F$ _( s( o
u(i)=b(i)/L(i);
6 c0 {, U# m\" A) N! H' m4 K, n
end7 s) {; h' k' f
L(N)=a(N)-c(N)*u(N-1);( e2 R4 {9 @3 ]0 L3 p& W
z(1)=d(1)/L(1);- U0 i& B: [8 h- T& N, p7 g
for i=2:N- k- [% W\" u, Z1 N
z(i)=(d(i)-c(i)*z(i-1))/L(i);
4 l4 d\" l' f/ K( e; ~; B5 y
end7 L+ f! n0 Y% }4 A
y(N)=z(N);1 g& `5 s: R1 H+ f& Y9 {9 L8 F( V
for i=N-1:-1:1* Q3 V) H u$ ]
y(i)=z(i)-u(i)*y(i+1);0 l S\" x+ P9 H6 c/ ^5 y$ h
end
! S7 B\" X- Y9 J' J
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%0 e9 H. J! n: _! I9 G/ }/ S# M! S
Y=[af,y,bt];5 `# y- s6 o8 [6 @- |/ j
for i=1:N+2
, O2 f, l8 H4 V
x=a0+(i-1)*h;
. X) N( l9 @/ `9 P5 C) n3 H+ y$ 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;: {- ]9 i, Y# w
end
: S4 }7 O# x4 ?; |
disp('下面两列分别是数值解和近似解');
% K' k; z& y+ D/ e: u* j
re=[Y' zj']