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

[复制链接]

字体大小: 正常放大

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');
* M\" V- e) S8 _! z
q=inline('2/x^2');
5 c! K. E: Q* w* a8 z% K
r=inline('sin(log10(x)/log10(exp(1)))/x^2'); k# Y7 U8 m7 E5 E9 T
N=9;+ F* s* B$ X, R' S: h `
a0=1;b0=2;
3 o5 c: |8 i* T2 H. q8 \
af=1;bt=2;\" t ]5 X+ g# d8 B5 g3 ]# U
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @$ ?! U- ~\" D
h=(b0-a0)/(N+1);) N9 m& ]. n& B( p7 K+ {! B
x=a0+h;
1 s( q: [2 I% N* m
a(1)=2+h*h*q(x); K6 h6 Z' E2 A/ x |
b(1)=-1+(h/2)*p(x);
# Z9 e8 T! G+ O: R1 S7 D# z7 @% {4 A
d(1)=-h*h*r(x)+(1+(h/2)*p(x))*af; U& p1 {; O9 C
for i=2:N-1
8 J3 @2 q% E5 w/ D5 Q; D\" c
x=a0+i*h;
1 c3 t3 x9 e- M0 }- [# ]. W
a(i)=2+h*h*q(x);
: I' C, H6 m) S* z/ y- d2 r
b(i)=-1+(h/2)*p(x);
* k6 t: c$ v1 E7 |* t
c(i)=-1-(h/2)*p(x);
6 n% G+ s\" a. P3 [, [: I
d(i)=-h*h*r(x);
& r# e) i. [. G7 }
end6 \- Q! j$ b( C y8 x
x=b0-h;) h$ ~& ]$ L$ H) _, R5 _9 U6 `7 i. B
a(N)=2+h*h*q(x);
2 y% Q9 a1 @/ _9 T2 e' y8 z
c(N)=-1-(h/2)*p(x);
) c6 D0 R7 d5 B, N8 |4 Y\" Y# y
d(N)=-h*h*r(x)+(1-(h/2)*p(x))*bt;
\" h8 e* e4 K8 B) T
%%%%%%%%%追赶法%%%%%%%%%%%%%%%%%%
; L7 t) C+ c* W, }( z
%y=trisys(c,a,b,d)5 k7 _3 L0 p' q! V% R3 ~
L(1)=a(1); K2 i\" e' F, F$ A
u(1)=b(1)/a(1);
) S, ^4 ]1 f% l4 [2 a
for i=2:N-19 M+ Z C) m$ Z( W
L(i)=a(i)-c(i)*u(i-1);, l2 Z% Q7 O1 X* H
u(i)=b(i)/L(i);, P5 _* z0 _3 p. p4 [1 y0 _
end9 K3 n' N* I8 K$ s3 L
L(N)=a(N)-c(N)*u(N-1);4 ^2 Z) t& R\" [, |- N
z(1)=d(1)/L(1);7 |- e. I( E1 |2 C6 w$ |
for i=2:N
d\" B+ n/ o0 x2 Z s
z(i)=(d(i)-c(i)*z(i-1))/L(i);
- P1 E) d1 Z& k: \, @
end
5 B$ g6 V4 n! @ V) ~+ q' }* q0 F
y(N)=z(N);
1 m% n0 j |7 ^* D3 l
for i=N-1:-1:1
( b) R. H% W+ R0 V, i
y(i)=z(i)-u(i)*y(i+1);
! i. P, u: Y; [7 C
end. u: |( @3 x$ c3 P& i
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%! U5 @9 N, H, s\" p* w
Y=[af,y,bt];
0 G) ] ]! x, h) f+ l2 Z
for i=1:N+2
! D$ I' ?0 E; m5 ^* N1 h8 _+ J\" D
x=a0+(i-1)*h;+ h3 n) A; V, R
zj(i)=1.1392070132*x-0.03920701320/x^2-3*sin(log10(x)/log10(exp(1)))/10-cos(log10(x)/log10(exp(1)))/10;8 _* `, W; L6 Z' j$ N
end9 a2 a( ?1 d6 [( Y' u) f
disp('下面两列分别是数值解和近似解');# b\" d& Z6 ?) I% b- x5 J
re=[Y' zj']