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

[复制链接]

字体大小: 正常放大

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');
, Q5 t$ Q# \1 g7 P' \- w; ^# b( v
q=inline('2/x^2');3 L4 N; ^7 d: z4 V2 G9 j
r=inline('sin(log10(x)/log10(exp(1)))/x^2');% Z1 [+ D& v! k6 {) S
N=9;7 I- e& s; T7 l9 O) O& z
a0=1;b0=2;
. e\" r: V7 D! [
af=1;bt=2;) M8 ]% g$ f5 m0 N: }. g# l
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
! _4 q/ P3 B: s( ^+ q) m+ g( I, Q
h=(b0-a0)/(N+1);7 n! z4 W. Z) O
x=a0+h;' u! W4 L* {* i/ L A+ x
a(1)=2+h*h*q(x);* M( I( ~, x2 n+ z\" R& m0 E
b(1)=-1+(h/2)*p(x);2 w5 O5 _9 T% U2 y a: ~
d(1)=-h*h*r(x)+(1+(h/2)*p(x))*af;
8 r1 Z! P: e6 u- P& x* P2 D
for i=2:N-1
# T3 B, i, e7 ?7 y5 B- [
x=a0+i*h;
1 ^6 h6 u: w% p) D* r3 l2 _3 \( d
a(i)=2+h*h*q(x);
* F5 K0 {; L0 ^6 L5 O
b(i)=-1+(h/2)*p(x);
\" H# \- Y! i' h9 L, `\" c' P* x9 R9 W
c(i)=-1-(h/2)*p(x);
' u. ^, K) Q7 I; R! m7 b
d(i)=-h*h*r(x);% _# v3 N% l6 r }$ x
end( X. ]: E4 k/ C\" r' E h! P
x=b0-h;' @* H3 \; U0 s5 H- j3 V+ i; L( Q
a(N)=2+h*h*q(x);
. ?. C. r0 f2 m0 N1 q\" x, n
c(N)=-1-(h/2)*p(x);+ j) Y4 M5 t5 O; _. g
d(N)=-h*h*r(x)+(1-(h/2)*p(x))*bt;2 O8 X. @ I/ x' b5 `6 u4 K
%%%%%%%%%追赶法%%%%%%%%%%%%%%%%%%
3 n b! d3 n\" v
%y=trisys(c,a,b,d)
5 S8 O- N\" W( x8 U P6 o\" V
L(1)=a(1);
# R3 o\" y+ j. w, B: t. c- H
u(1)=b(1)/a(1);2 z2 U% L' L4 L7 M
for i=2:N-11 _\" C2 M1 W ]1 ?9 T7 A* P
L(i)=a(i)-c(i)*u(i-1);
4 L\" A7 {& N5 A! e( k/ x) c5 `
u(i)=b(i)/L(i);
, X; N% s2 t; D) e9 e3 i
end1 v8 i1 _7 d5 m3 r$ i2 Y. { }
L(N)=a(N)-c(N)*u(N-1);6 _5 I; N$ Y, l( @
z(1)=d(1)/L(1);
8 ?\" G( Z9 {. n9 H
for i=2:N6 K8 i! s. q: Z8 w\" Y- U
z(i)=(d(i)-c(i)*z(i-1))/L(i);
4 |! q7 l2 e1 F: v0 }
end& f+ [) k0 P8 ]/ p' g
y(N)=z(N);
5 w( A+ ~/ _; J- S' m
for i=N-1:-1:1
8 A* ~' p$ [* ~, ?
y(i)=z(i)-u(i)*y(i+1);
* @ Q\" E. P% j+ P Y2 v) M
end8 T# p\" m/ [2 q8 o
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# Q\" @8 i4 M& u, _' O
Y=[af,y,bt];! C7 B: @! S% D
for i=1:N+2
/ @% G; ^8 h. c3 L/ g4 c
x=a0+(i-1)*h;
' F& M5 z, N7 ?6 p2 }: D9 J
zj(i)=1.1392070132*x-0.03920701320/x^2-3*sin(log10(x)/log10(exp(1)))/10-cos(log10(x)/log10(exp(1)))/10;6 X% f4 W* G! N$ m, T
end
( W( Z4 E4 ]9 ^. `\" M
disp('下面两列分别是数值解和近似解');
3 \6 e4 X2 s- N5 e\" o/ P G( K& p
re=[Y' zj']