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

[复制链接]

字体大小: 正常放大

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');
4 t) P% h( w( P3 `
q=inline('2/x^2');
0 Q\" V7 B o+ d
r=inline('sin(log10(x)/log10(exp(1)))/x^2');
r* f8 D) X* K' U7 T A
N=9;
8 x; C6 _- Y- y+ {/ w
a0=1;b0=2;0 u* a. Z0 S% @! f1 m
af=1;bt=2;3 g: \ N: M, n* S; X( V
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%9 Q2 p\" R5 V& l1 a7 |, m\" D6 _
h=(b0-a0)/(N+1);2 S I: B2 z/ O! X B
x=a0+h;* t ?1 U% N7 {; m, t% T
a(1)=2+h*h*q(x);
n+ V$ V; ?* o# j
b(1)=-1+(h/2)*p(x); r- V- n& B6 B/ p
d(1)=-h*h*r(x)+(1+(h/2)*p(x))*af;% `, j. D* N% x9 B, z
for i=2:N-1
/ L7 r% H* \; U/ p
x=a0+i*h;+ n\" |# U; `' q1 j
a(i)=2+h*h*q(x);6 [6 @3 S% Q2 Z\" G$ k5 @! U\" _) ]
b(i)=-1+(h/2)*p(x);
9 \) @) w; b( n
c(i)=-1-(h/2)*p(x);
( ~: z* a) `& O% M! S0 }/ B
d(i)=-h*h*r(x);+ S. \2 X6 k/ S$ K- ^5 i* l9 O
end, Q' A! Z: w! Z1 V
x=b0-h;
* L3 y! l) h6 W! y7 x
a(N)=2+h*h*q(x);2 t E5 {9 z, v% F! ?. a/ j
c(N)=-1-(h/2)*p(x);; j; [/ L9 [& l+ c
d(N)=-h*h*r(x)+(1-(h/2)*p(x))*bt;) }/ e2 j4 O8 e' D\" L, t ^9 p
%%%%%%%%%追赶法%%%%%%%%%%%%%%%%%%
\" p: j& x( c% R! m+ I
%y=trisys(c,a,b,d)/ h) s$ w) ^! T2 u! q
L(1)=a(1);3 c- j9 i( Z, C! b3 N9 F. W$ t
u(1)=b(1)/a(1);
) E3 {( s6 D- C$ ^+ ]
for i=2:N-1
3 q$ @4 D) \: ^+ ^+ m q: `
L(i)=a(i)-c(i)*u(i-1);9 P& M1 q/ s! d! y9 H
u(i)=b(i)/L(i);/ I0 H. D& U' w7 B( E
end
( Y; A- L$ Y, g6 P. B, ^& A
L(N)=a(N)-c(N)*u(N-1);# C, {7 I3 d8 s
z(1)=d(1)/L(1);+ v& P/ l# x% N8 z
for i=2:N/ P& m; t0 H/ [7 `\" D
z(i)=(d(i)-c(i)*z(i-1))/L(i);( N& ^( j/ C# B( {2 R' i
end
0 B8 {( ?% n) p, B% l4 o3 l6 q
y(N)=z(N);: L8 {, q7 E ^) \\" D) N1 c m
for i=N-1:-1:14 E0 i: Y, o' F. F' a
y(i)=z(i)-u(i)*y(i+1);9 u; i) _: \1 v. K1 d# A2 ^
end/ w% y( S% q; w- d\" O2 B. \; v
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5 e: a2 ]$ s& q1 ?1 x& o) C
Y=[af,y,bt];$ v: l* L- p; W\" n4 i! r; R
for i=1:N+2
; U, N2 v\" Z+ q) n# A9 S
x=a0+(i-1)*h;
8 Z, a7 _) ?\" v- k* Q/ W$ o
zj(i)=1.1392070132*x-0.03920701320/x^2-3*sin(log10(x)/log10(exp(1)))/10-cos(log10(x)/log10(exp(1)))/10;
5 b\" j7 D: o* M5 I7 N, ~0 h3 ]
end5 K; k* v: Z9 M+ o4 Q3 ?
disp('下面两列分别是数值解和近似解');5 L% z& e1 N' `: k' u1 u3 j' w
re=[Y' zj']