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

[复制链接]

字体大小: 正常放大

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');' Q1 J; `+ m\" J2 o: ^
q=inline('2/x^2');/ G\" P3 W' B* c
r=inline('sin(log10(x)/log10(exp(1)))/x^2');* Y! b. y/ @3 J3 {7 t3 D
N=9;+ q+ M5 F; A, E$ \8 `
a0=1;b0=2;
5 R: o4 |! A2 q- `% S5 c
af=1;bt=2;
G' p \# i4 P7 m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%: m& X! s# I) O3 ?
h=(b0-a0)/(N+1);
; s4 j7 y( [3 ^* ^; c3 I% K
x=a0+h;& [; k- w1 W4 j9 Z1 p1 i
a(1)=2+h*h*q(x);
8 j0 z& S/ Z& S6 y6 l; g' m' Z
b(1)=-1+(h/2)*p(x);* A$ V! {3 _! `- T, l, Y$ z
d(1)=-h*h*r(x)+(1+(h/2)*p(x))*af;- i/ _, u3 x& |4 d0 J/ U% d
for i=2:N-1% Q# n\" w$ i( p3 X {7 c9 M
x=a0+i*h;
: K1 d+ k! F' B+ ^/ _# W
a(i)=2+h*h*q(x);
& j: a! f. F2 m) e4 q7 ]
b(i)=-1+(h/2)*p(x);
4 `, w0 T7 o+ }4 L7 G5 @' r3 E) R X
c(i)=-1-(h/2)*p(x);
# `! |7 E: K& L& h5 E
d(i)=-h*h*r(x);3 `4 N8 z0 x& A
end0 E, E3 D: g% }- C
x=b0-h;. ~- T' Q. O% i1 J, y6 C9 a
a(N)=2+h*h*q(x);; O1 P! @% E) p& ^6 B
c(N)=-1-(h/2)*p(x);
7 i' s9 D3 i# D8 t* V* {, J\" e
d(N)=-h*h*r(x)+(1-(h/2)*p(x))*bt;
9 [& L/ r6 f1 V( y# B( Q* P- I/ j
%%%%%%%%%追赶法%%%%%%%%%%%%%%%%%%- [2 |0 k4 f\" B; D' Z9 a6 M
%y=trisys(c,a,b,d)
# Q9 u- Z4 y5 h+ _2 |( T4 n/ V
L(1)=a(1);
' u% L p$ n i9 \1 ]: S
u(1)=b(1)/a(1);
I5 A% s# r. S* Q( \* l. o7 [
for i=2:N-1
$ c, h: i5 i! G; G6 w
L(i)=a(i)-c(i)*u(i-1);
* f9 E- ~+ i% I9 s5 A. U4 z G
u(i)=b(i)/L(i);# i' T4 n: ]) R! c% G) a+ ~
end) g2 p1 g: g2 m+ t: C8 e
L(N)=a(N)-c(N)*u(N-1);) ~: J7 _, Y( L9 A0 r8 _. {
z(1)=d(1)/L(1);
+ y; p; N% _3 u8 w0 |4 s
for i=2:N
9 V# I6 {6 p6 B: D
z(i)=(d(i)-c(i)*z(i-1))/L(i);# t5 Q$ u( G, k% c6 h- Q9 B
end+ ^& F) h3 n0 M% w
y(N)=z(N);0 W* b/ \; g; l
for i=N-1:-1:1, z% ]5 B+ t `; F
y(i)=z(i)-u(i)*y(i+1);( h. h; h0 o# X
end
+ Z& j3 d* V& L) |( r% X
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%# R( n- j. g% D K( I# F' T- K5 g* L
Y=[af,y,bt];
% k$ \; K( t5 j
for i=1:N+2
/ N\" m* x# F$ V) n
x=a0+(i-1)*h;
6 A; e9 g* s& l a6 r) x
zj(i)=1.1392070132*x-0.03920701320/x^2-3*sin(log10(x)/log10(exp(1)))/10-cos(log10(x)/log10(exp(1)))/10;- ^7 G, w& N4 V3 U3 u
end
* \6 e a& `4 J
disp('下面两列分别是数值解和近似解');- l) m) [; v9 s* S7 b( {% i7 d& c* T: i
re=[Y' zj']