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

[复制链接]

字体大小: 正常放大

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');) G/ [' A1 Z1 `\" j9 h# j7 r
q=inline('2/x^2');! P% g- ?( u\" N+ I/ q
r=inline('sin(log10(x)/log10(exp(1)))/x^2');
, A d( E2 w$ q# m6 |; T
N=9;
5 p! [1 F* {5 ]2 N; f) o
a0=1;b0=2;
7 Y/ C' e2 t( \$ x: B+ T o/ Y3 I
af=1;bt=2;+ K4 ?, ]2 I$ m. T8 n; G; {: j
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%0 w' d3 S( X5 H; ~
h=(b0-a0)/(N+1);
% j, N( z# ~: g
x=a0+h;
; C4 K% F* P3 }! p: T: |( B
a(1)=2+h*h*q(x);
5 I7 m3 l! }1 S- ]. u
b(1)=-1+(h/2)*p(x);
! W& l, @& }6 p) ?2 _
d(1)=-h*h*r(x)+(1+(h/2)*p(x))*af;5 S% K0 r5 o6 E6 q+ B4 u9 @$ f, b$ u\" m
for i=2:N-1
x=a0+i*h;
# L) @( a% k; `
a(i)=2+h*h*q(x);% W* |4 W* Z\" f# p
b(i)=-1+(h/2)*p(x);; O3 [0 a$ q+ z4 F( H& h
c(i)=-1-(h/2)*p(x);
\" {- R\" K; K/ a2 I A
d(i)=-h*h*r(x);; o2 W* W) U! v: V4 t% I: T
end* q; x' w/ m) O! b5 z\" a- d
x=b0-h;! \) r0 N' Z) F; e5 c
a(N)=2+h*h*q(x);$ e' l( P4 ?8 f- K5 _) G3 u
c(N)=-1-(h/2)*p(x);: _+ |& a! ~2 A9 ]5 S7 a8 a
d(N)=-h*h*r(x)+(1-(h/2)*p(x))*bt;3 t- B w& c0 S) i& r! ^4 _8 u8 ?
%%%%%%%%%追赶法%%%%%%%%%%%%%%%%%%
\" s# X' O6 F) ?1 }# Q$ h
%y=trisys(c,a,b,d): X& g* b! X+ k. ~: C# A& \& I+ [
L(1)=a(1);
9 K+ _, @- U3 W4 U! T# s z
u(1)=b(1)/a(1);9 I2 b1 D% y% h/ x+ R
for i=2:N-1
+ Z- a( Y0 u4 b; ?0 A1 x, K
L(i)=a(i)-c(i)*u(i-1);
0 o& q- p: B, O, \, ^
u(i)=b(i)/L(i); y( O- {, t2 N
end$ V! |1 I% A. B+ G
L(N)=a(N)-c(N)*u(N-1);+ n0 x) b( W- v* o& o
z(1)=d(1)/L(1);
7 _9 {9 h3 X0 S8 A
for i=2:N: I7 t% D D9 M) n- o! h6 ]9 v
z(i)=(d(i)-c(i)*z(i-1))/L(i);
% z) g1 s- K! p9 o
end! i6 C4 F9 W# F1 r9 a/ P+ H
y(N)=z(N); x8 s: I5 E1 c6 C/ d0 c: A
for i=N-1:-1:1
, B( j5 q) j6 p3 x
y(i)=z(i)-u(i)*y(i+1);
; D7 Q8 v5 w$ X# g) g9 }
end
1 t8 V3 i0 }$ p! r8 ?. A
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% I# v! B2 ?, T. {5 b2 a
Y=[af,y,bt];4 b- V5 [+ z6 h; T
for i=1:N+2. G; \& H0 C8 L9 r
x=a0+(i-1)*h;
- o# }( |7 e. e' J% 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;7 w* l* l8 a6 i# a
end: ~- ]8 R+ o' t% K. V$ h- e
disp('下面两列分别是数值解和近似解');) T+ S: E1 m, N+ h
re=[Y' zj']