有限差分法和托马斯算法（或追赶法）对一个二阶线性边值问题进行数值求解 - 数学建模社区-数学中国

p=inline('-2/x');
4 C8 ?/ q& g+ H: k! r
q=inline('2/x^2');
) i2 p& @) e! F8 U) Z
r=inline('sin(log10(x)/log10(exp(1)))/x^2');* ]0 Y9 p. t/ M1 W3 X- L% n) |1 x
N=9;
$ Q+ H9 o2 A* h2 {& x
a0=1;b0=2;: v3 v! Z7 K% ^4 N! i
af=1;bt=2;
X) {: ^5 r0 c" P7 \5 T: I
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%# p/ e) m0 _. h+ x
h=(b0-a0)/(N+1);+ w7 g/ b' U5 p
x=a0+h;
. P8 ?0 V" ?4 w$ n9 q2 r6 j( M6 P
a(1)=2+h*h*q(x);
3 w; n7 J, ]+ ~% k$ s1 @4 b
b(1)=-1+(h/2)*p(x);
5 b# P9 [! \' j! B/ b+ m
d(1)=-h*h*r(x)+(1+(h/2)*p(x))*af;
# t* R f! n7 y' F% ~8 m
for i=2:N-1
" V+ U- n3 u' \6 V [
x=a0+i*h;
- j. o. j6 S& A2 ?# f6 O
a(i)=2+h*h*q(x);: t9 V7 r4 N# m( }$ a0 J* A
b(i)=-1+(h/2)*p(x);2 q. K5 i* F6 R/ h( @" j# s
c(i)=-1-(h/2)*p(x);. ^6 K8 ~; }" F- @
d(i)=-h*h*r(x);
' }' c7 f$ R8 H' m0 ^
end- h9 u* M) P0 \ o4 f( ^/ M
x=b0-h;
" H! U4 |$ O8 @, l. E" q
a(N)=2+h*h*q(x);( b; G( t1 ]3 D- @* y5 f6 u
c(N)=-1-(h/2)*p(x);. q7 J @ X# T. F
d(N)=-h*h*r(x)+(1-(h/2)*p(x))*bt;' W5 |/ G. d3 W5 _
%%%%%%%%%追赶法%%%%%%%%%%%%%%%%%%
+ {" c/ h. I( ~: ~; b$ X2 R
%y=trisys(c,a,b,d)' H- R8 ?$ ]9 B a' k, z( Y& o
L(1)=a(1);
5 t2 C+ k; `$ U s; {
u(1)=b(1)/a(1);! _$ a* v# o0 f* `
for i=2:N-17 B' q; v1 h9 Z [4 Y
L(i)=a(i)-c(i)*u(i-1);0 \# n- B, U o/ N! z' j3 {& U$ ~
u(i)=b(i)/L(i);
7 U- `+ l/ i2 g
end
7 I2 O- s/ T5 Q/ ]3 A) v- a v
L(N)=a(N)-c(N)*u(N-1);
' A( d) t; ]( g C$ N
z(1)=d(1)/L(1);
1 R, B$ H3 O. c$ Y
for i=2:N
8 M0 B% }, B9 [2 M1 C
z(i)=(d(i)-c(i)*z(i-1))/L(i);+ ^% h8 @) F0 p7 M) a" a
end
/ U: n$ O# |+ a) L! a
y(N)=z(N);
! b) S2 U" V- C! q+ a! _1 U
for i=N-1:-1:19 D$ e! c% J) {
y(i)=z(i)-u(i)*y(i+1);$ ?. s: b' \) y
end
% U1 a+ d9 k6 _! d" F" w3 [
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 x7 _* j {1 G, h% f
Y=[af,y,bt];
; q6 N1 \5 m5 u ^. I3 {5 }
for i=1:N+2
: }4 _: F3 [+ v
x=a0+(i-1)*h;$ B; t3 w. u5 I
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 Z% [0 |. \; z+ r+ y: ~ X
end
# v J+ [3 \. U/ y( k# C# J& ]
disp('下面两列分别是数值解和近似解');
% G' ~/ ^; p& U- z. [# E" W
re=[Y' zj']

复制代码