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

[复制链接]

字体大小: 正常放大

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');# J, p: p: ^- v; a+ S
q=inline('2/x^2');
+ N* H0 t) r% Z( Z- m+ v1 c& ~
r=inline('sin(log10(x)/log10(exp(1)))/x^2');. I2 `7 |4 ^) M/ {- d
N=9;2 T7 U& O# q' d1 H
a0=1;b0=2;% y( B0 I- {' v W, j3 v
af=1;bt=2;1 r3 i ]5 u4 T3 _# a
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%1 `' I- I r' Y; e4 O
h=(b0-a0)/(N+1);
e$ i' _2 J: r5 m X
x=a0+h;. [4 B( R\" _\" l$ K, n
a(1)=2+h*h*q(x);
8 f8 e- `0 x\" b: a\" } i
b(1)=-1+(h/2)*p(x);
_5 q) \7 b: x& r, Z1 t8 }; h: m
d(1)=-h*h*r(x)+(1+(h/2)*p(x))*af;
) s- X; q\" N; `2 x/ H2 v
for i=2:N-19 B+ A, z- O: s\" ]
x=a0+i*h;1 J- W% x% ]: ~* Y( _
a(i)=2+h*h*q(x);# m+ s9 X$ J( g8 A1 J v8 E
b(i)=-1+(h/2)*p(x);
( ~, P+ q8 o/ M. M7 U6 R# @
c(i)=-1-(h/2)*p(x);
, R9 J* L1 P8 O. G$ h2 {+ h M9 m
d(i)=-h*h*r(x);
& E6 V) Y. i; N/ [0 o
end$ D# \\" K# a q) R( w' |
x=b0-h;8 z* R( g/ ?\" f( H
a(N)=2+h*h*q(x);
n% W T B, A% e9 L0 b
c(N)=-1-(h/2)*p(x);! G\" Y' G# U. o% r
d(N)=-h*h*r(x)+(1-(h/2)*p(x))*bt;
1 N: f F- A7 x5 F1 _& {& O. \
%%%%%%%%%追赶法%%%%%%%%%%%%%%%%%%; a& Z- |$ I9 m. ^' q, q4 P3 s
%y=trisys(c,a,b,d)
0 c: A3 L% h P4 R/ `2 b) ^
L(1)=a(1);
, r0 a+ p3 N* o% ~9 m& K; P/ |5 E
u(1)=b(1)/a(1);- j5 b2 S$ R/ z5 b( q c1 R& J
for i=2:N-1
) O* B% G( a+ B8 ^9 {
L(i)=a(i)-c(i)*u(i-1); k0 @+ c! R/ Y& l
u(i)=b(i)/L(i);4 k/ l& q/ N1 X
end
2 H, P, r; o$ u8 U
L(N)=a(N)-c(N)*u(N-1);
1 U' Z, r4 T& Q! ]$ K0 c
z(1)=d(1)/L(1);- p4 n/ r' t& U\" ^
for i=2:N
0 c$ U7 B4 P; d\" a) W
z(i)=(d(i)-c(i)*z(i-1))/L(i);! z/ m2 h. X/ C+ ^3 ^
end
# ]( _3 K8 L* Y/ K% Q
y(N)=z(N);
\" ]: @. X* j; g: ~ j\" ^
for i=N-1:-1:17 ~5 l6 L9 w4 t/ Z
y(i)=z(i)-u(i)*y(i+1);( `1 ~: C* s4 @, x
end
$ p9 m0 i) M* F' C: W2 W
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%3 `1 u9 K9 ]/ j! Q
Y=[af,y,bt];6 r, u% Q, n! C+ q9 P
for i=1:N+2
2 B7 t6 h\" a0 D
x=a0+(i-1)*h;, _/ @- U. z( C8 F+ R% k
zj(i)=1.1392070132*x-0.03920701320/x^2-3*sin(log10(x)/log10(exp(1)))/10-cos(log10(x)/log10(exp(1)))/10;! S/ T6 B! ?2 n# c( M
end9 Y5 n$ C3 |) \# H4 u) p
disp('下面两列分别是数值解和近似解');
+ @' M/ H$ e7 P
re=[Y' zj']