科学与工程计算编程作业1

科学与工程计算作业

考虑如下方程：

220002201

1,(0,1)0

10(()),()(10)10t t x x x x u u x t x

u u

v t u q v t q dt x

u

u x =====??=∈??=?=-==-??=?? 要求如下：

1：两步格式、紧格式、三层隐格式、跳点格式进行计算。（宇航、机电、车辆） 2：最简显、隐格式、Crank-Nicolson 格式、DFF 格式进行计算。（其他学院）

解：

1、矩形网格剖分区域

根据题意要求：将变量x 的计算区间[0,1]进行区间分割如下：

取空间步长h =110

, 时间任意的矩形网格剖分区域, 用节点表示坐标点(j,k)表示坐标点1(,)(,k )j 0,1,...,,k 0,1,...,j k T x t jh h ττ===,

0101M x x x =???= 1 , , 1 , 0

j J k h x jh x Mh M t k n τ=====≥

以后沿用记号为()k ,t k

j j u u x =

作为经常使用的标记。

2、建立差分格式

(1)最显简格式： 变形有:1112(12)k k k k

j j j j u ru r u ru r h τ+-+??

=+-+=

??

?

(2)最显隐格式：1112

11（）（2）k k k k k

j j j j j u u u u u h

τ-+--=-+

变形有：1

11(2)k k k k k j j j j j u r u u u u --+--+=

(3) Crank-Nicolson 格式:

1111

11112

1（[2）（2）]2k k k k k k k k i i i i i i i i u u u u u u u u h

τ+++++-+--=-++-+ 变形有：11111111（[2）（2）]2

k k k k k k k k i i i i i i i i r

u u u u u u u u +++++-+-=+

-++-+

(4) Du-Fort 和Frankel （DFF ）格式：

11111

1

2

()02τ

+-+-+---++-

=k k k k k k j

j

j j

j

j u

u

u

u

u u

h

变形有：

11111

12(())

+-+-+-=+-++k k k k k k

j

j

j j

j

j u

u

r u

u

u u

3、边界处置离散问题，该题属于第二类或第三类初边值问题

初值离散0

00j u j M =≤≤

边界的离散可以按下面的方法进行：

(1)在点(0,)和点(1,)分别用向前和向后差商近似，有：

k k t t u x ??10

1()()

k k

k

k k

k

k M M M k

u u u t h

u u u t h

αμβγ-?-=+???-?=+??

22(,)(,)0,01i k i k u u

x t x t x t x ??-=<

1211（）（2）

k k k k k i i i i i u u u u u h

τ++--=-+

与内点查分格式一起，即可求解，但边界条件的离散是一阶的，不适用于二阶格式。

（2）在边界节点上的一阶导数可以构造一阶中心差分格式，即

这种离散是二阶的，E=O(h 2),但式中含有区域之外的点 ,需要借助内点差分格式消去11,k k

M u u -+

例如：内点计算用古典显示格式 分别取j=0和j=M 得：

与前面的离散格式联立，得到如下公式：

这样边界节点的值可以通过内部节点的值表示出来，而且其误差精度是二阶的。和内部节点所得到的方程一起联立求解。

根据题意所给的要求，可以到： 110,10,1,0.1,4()0,()()

M k a h r t u t v t αβγ=-====

== 则上式化简为：10011111()22012

++-?=+-????=+??k k k k k k k M M M u u u v t u u u 4、求解题意中给出的v(t),有题提供的

0022010(()),()(10)==?=-==-??t x x u v t u q v t q dt x 可以得出 110110()2()2k k k k k k k M M k u u u t h u u u t h αμβγ-+-?-=+???-?=+??1

1(,),(,)k M M x t x t -+111220k k k k k j j j j j u u u u u a h

τ++---+-=100101111(2)(2)k k k k k k k k k k M M M M M u u ar u u u u u ar u u u +-++-=+-+=+-+10011112(1)22()12(1)22()k k k k k k k M

M M M u ar ah u aru arhu t u ar h u aru arh t βγ++-???=-++-??????=--++????

00()(10(10(())))==--?t x v t v t u dt

有题意知(0)0(0,0)0，==v u ,上面的式子可写出如下微分方程

10(1()(0,))=-+dv v t u t dt

根据以前所学的常微分方程数值解法，可以利用四阶龙格—库塔方法(R-K 方法)

1123411

22

3431

(22)

6(,)

(,)

22(,)

22(,)

0,1,2,...,10

+=++++==++=++=++=n n n n n n n n n n v v k k k k k hf t v k h k hf t v k

h k hf t v k hf t h v k n

则知道边界条件的离散值：

1001110109

1123411

22

34300

011

()

2201

21

(22)

6(,)

(,)

22(,)

22(,)

0,1,2,...

τττττττ+++?=+-????=+??=++++==++=++=++===k k k k

k k k

n n n n n n n n n n u u u v t u u u v v k k k k k f t v k k f t v k k f t v k f t v k n v u 5、几种格式的差分方程组：

（1）最简隐格式差分方程组

1110100111010911()0,1,2,,10 k 0,1,2,2400,1,2,,1011()k 0,1,2,2201k 0,1,2,2

++-++?=++=???=?????==??????=+-=??????=+=?????k k k k j j j j j k k k k k k k u u u u j u j u u u v t u u u

具体的Matlab 的Expscheme.m 文件的代码如下

f=@(t,v,u)(10*(1-v+u));%求v(t)函数的句柄。

h=1/10;

T = 0.5;

tau= 0.005;

M=1/h;

N=T/tau;

x = [0:M]'*h;

t = [0:N]*tau;

r = tau/h/h ；

r1 = 1 - 2*r;

u=zeros(M+1,1);

V(1)=0;

for k = 1:N

k1=tau*feval(f,t(k),V(k),u(1,k));

k2=tau*feval(f,t(k)+tau/2,V(k)+k1/2,u(1,k));

k3=tau*feval(f,t(k)+tau/2,V(k)+k2/2,u(1,k));

k4=tau*feval(f,t(k)+tau,V(k)+k3,u(1,k));

V(k+1)=V(k)+(k1+2*k2+2*k3+k4)/6;

u(1,k+1)=u(1,k)+0.5*u(2,k)-1/20*V(k);

for i = 2:M %?ó?ú2??μ

u(i,k+1) = r*(u(i + 1,k) + u(i-1,k)) + r1*u(i,k);

end

u(M+1,k+1)=u(M+1,k)+0.5*u(M,k);

end

surf(t,x,u);

xlabel('t')

ylabel('x')

zlabel('u')

title('最简显格式误差,时间步长为1/200,空间步长为1/10’)

当10.010.510

T h τ===

（2）最简隐格式差分方程组

1111010011101091(2)0,1,2,,10 k 0,1,2,400,1,2,,1011()k 0,1,2,2201k 0,1,2,2

+-+-++?=+-+=???=?????==??????=+-=??????=+=?????k k k k k j j j j j j k k k k k k k u u u u u j u j u u u v t u u u 具体的Matlab 的Impscheme.m 文件的代码如下

f=@(t,v,u)(10*(1-v+u));%?óv(t)oˉêyμ???±ú￡???±?á??at

h=1/10;

T = 0.5;

tau= 0.005;

M=1/h;

N=T/tau;

x = [0:M]'*h;

t = [0:N]*tau;

r = tau/h/h, r1 = 1 - 2*r;

u=zeros(M+1,1);

V(1)=0;

for k = 1:N

k1=tau*feval(f,t(k),V(k),u(1,k));

k2=tau*feval(f,t(k)+tau/2,V(k)+k1/2,u(1,k));

k3=tau*feval(f,t(k)+tau/2,V(k)+k2/2,u(1,k));

k4=tau*feval(f,t(k)+tau,V(k)+k3,u(1,k));

V(k+1)=V(k)+(k1+2*k2+2*k3+k4)/6;

u(1,k+1)=u(1,k)+0.5*u(2,k)-1/20*V(k);%?ó×ó±??μ for i = 2:M %?ó?ú2??μ

u(i,k+1)=u(i,k)+1/4*(u(i-1,k) -2*u(i,k)+ u(i+1,k)); end

u(M+1,k+1)=u(M+1,k)+0.5*u(M,k); %?óóò±??μ

end

surf(t,x,u);

xlabel('t')

ylabel('x')

zlabel('u')

title('最简隐格式误差,时间步长为1/200,空间步长为1/10')

当10.010.510

T h τ===

（3）Crank-Nicolson 格式差分方程组

11111111010011101091[22]0,1,2,,10 k 0,1,2,800,1,2,,1011()k 0,1,2,2201k 0,1,2,2（）（）+++++-+-++?=+-++-+=???=???

??==??????=+-=??????=+=?????k k k k k k k k j i i i i i i i j k k k k k k k u u u u u u u u j u j u u u v t u u u 具体的Matlab 的CrankNicolson .m 文件的代码如下

f=@(t,v,u)(10*(1-v+u));%?óv(t)oˉêyμ???±ú￡???±?á??at

h=1/10;

T = 0.5;

tau= 0.005;

M=1/h;

N=T/tau;

x = [0:M]'*h;

t = [0:N]*tau;

r = tau/h/h, r1 = 1 - 2*r;

u=zeros(M+1,1);

V(1)=0;

for k = 1:N

k1=tau*feval(f,t(k),V(k),u(1,k));

k2=tau*feval(f,t(k)+tau/2,V(k)+k1/2,u(1,k));

k3=tau*feval(f,t(k)+tau/2,V(k)+k2/2,u(1,k));

k4=tau*feval(f,t(k)+tau,V(k)+k3,u(1,k));

V(k+1)=V(k)+(k1+2*k2+2*k3+k4)/6;

u(1,k+1)=u(1,k)+0.5*u(2,k)-1/20*V(k);

for i = 2:M

u(i,k+1)=u(i,k)+0.125*r*((u(i -1,k) -2*u(i,k)+ u(i+1,k))+(u(i -1,k+1) -2*u(i,k+1)+ u(i+1,k+1)));

end

u(M+1,k+1)=u(M+1,k)+0.5*u(M,k);

end

surf(t,x,u);

xlabel('t')

ylabel('x')

zlabel('u')

title('CrankNicolson 格式误差,时间步长为1/200,空间步长为1/10')

当10.010.510

T h τ===

（3）DFF 格式差分方程组

111111010011101091(())0,1,2,,10 k 0,1,2,200,1,2,,1011()k 0,1,2,2201k 0,1,2,2+-+-+-++?=+-++=???=?????==??????=+-=??????=+=?????k k k k k k j j j j j j j k k k k k k k u u u u u u j u j u u u v t u u u

Matlab 对于的代码为：DFF.m

f=@(t,v,u)(10*(1-v+u));%?óv(t)oˉêyμ???±ú￡???±?á??at

h=1/10;

T = 0.5;

tau= 0.005;

M=1/h;

N=T/tau;

x = [0:M]'*h;

t = [0:N]*tau;

r = tau/h/h, r1 = 1 - 2*r;

u=zeros(M+1,1);

V(1)=0;

for k = 1:N

k1=tau*feval(f,t(k),V(k),u(1,k));

k2=tau*feval(f,t(k)+tau/2,V(k)+k1/2,u(1,k));

k3=tau*feval(f,t(k)+tau/2,V(k)+k2/2,u(1,k));

k4=tau*feval(f,t(k)+tau,V(k)+k3,u(1,k));

V(k+1)=V(k)+(k1+2*k2+2*k3+k4)/6;

u(1,k+1)=u(1,k)+0.5*u(2,k)-1/20*V(k);

for i = 2:M

u(i,k+2)=u(i,k+1)+0.5*(u(i+1,k+2)-(u(i,k+2)+u(i,k))+u(i-1,k+1)); end

u(M+1,k+1)=u(M+1,k)+0.5*u(M,k);

end

surf(t,x,u);

xlabel('t');

ylabel('x');

zlabel('u');

title('DDF 差分格式时间误差为1/200,空间误差1/10’);

0.25

t r=0.25

x

u

本文来源：https://www.bwwdw.com/article/ua6l.html