微分方程作业.docx

微分方程作业.docx

《微分方程作业.docx》由会员分享，可在线阅读，更多相关《微分方程作业.docx（12页珍藏版）》请在冰豆网上搜索。

微分方程作业

P10习题

1.用Euler法和改进的Euler法求u’=-5u（0≤t≤1）,u（0）=1的数值解，步长h=0.1,0.05;并比较两个算法的精度。

解：

functiondu=Euler_fun1（t,u）

du=-5*u;clear;

h=0.1;tend=1;N=1/h;t

（1）=0;u

（1）=1;

t=h.*（0:

N）;

forn=1:

N

u（n+1）=u（n）+h*Euler_fun1（t（n）,u（n））;

end

plot（t,u,'*'）;holdon

forn=1:

N

v

（1）=u（n）+h*Euler_fun1（t（n）,u（n））;

fork=1:

6

v（k+1）=u（n）+h/2*（Euler_fun1（t（n）,u（n））+Euler_fun1（t（n+1）,v（k）））;

end

u（n+1）=v（k+1）;

end

plot（t,u,'o'）;

sol=dsolve（'Du=-5*u','u（0）=1'）;

u_real=eval（sol）;

plot（t,u_real,'r'）;

将上述h换为0.05得：

由图像知道:

显然改进的Euler法要比Euler法精确度要高；

3.将u‘’=-u（0≤t≤1），u（0）=0,u’（0）=1化为一阶方程组，并用Euler法和改进的的Euler法求解，步长h=0.1,0.05；并比较两个算法的精度。

解：

functiondu=fun31（y）

du=y;

functiondy=fun32（u）

dy=-u;

clear;

h=0.1;tend=1;N=1/h;t

（1）=0;u

（1）=0;y

（1）=;

t=h.*（0:

N）;

forn=1:

N

u（n+1）=u（n）+h*y（n）;

y（n+1）=y（n）+h*（-u（n））;

end

plot（t,u,'*'）;holdon

forn=1:

N

v

（1）=u（n）+h*fun31（y（n））;

w

（1）=y（n）+h*fun32（u（n））;

fork=1:

6

v（k+1）=u（n）+h/2*（fun31（y（n））+fun31（...w（k）））;

w（k+1）=y（n）+h/2*（fun32（u（n））+fun32（...v（k）））;

end

u（n+1）=v（k+1）;

y（n+1）=w（k+1）;

end

plot（t,u,'o'）;

sol=dsolve（'D2u=-u','u（0）=0','Du（0）=1';

u_real=eval（sol）;

plot（t,u_real,'r'）;

将上述h换为0.05得：

由图像可以知道:

显然改进的Euler法要比Euler法精确度要高；

实习题

（二）

1.取步长

，分别用Euler法和改进的Euler法求下列初值问题的解，并与真解相比较.

（1）

真解

；

解：

functiondu=fun1（x,u）

du=u-2*x/u;

clear;

h=0.1;xend=1;N=1/h;x

（1）=0;u

（1）=1;

x=h.*（0:

N）;

%——Eluer法——%

forn=1:

N

u（n+1）=u（n）+h*fun1（x（n）,u（n））;

end

plot（x,u,'*'）;

holdon

%——改进的Eluer法——%

forn=1:

N

v

（1）=u（n）+h*fun1（x（n）,u（n））;

fork=1:

6

v（k+1）=u（n）+h/2*（fun1（x（n）,u（n））+fun1（x（n+1）,v（k）））;

end

u（n+1）=v（k+1）;

end

plot（x,u,'o'）;

holdon

%——真解——%

u_real=sqrt（1+2*x）;

plot（x,u_real,'r'）;

由图像可以知道:

显然改进的Euler法要比Euler法精确度要高；

（2）

真解

；

解：

functiondu=fun2（x,u）

du=（u/x）-x.^2/u.^2;

clear;

h=0.1;N=1/h;x=1:

h:

2;x

（1）=1;u

（1）=2;

forn=1:

N

u（n+1）=u（n）+h*fun2（x（n）,u（n））;

end

plot（x,u,'*'）;

holdon

forn=1:

N

v

（1）=u（n）+h*fun2（x（n）,u（n））;

fork=1:

6

v（k+1）=u（n）+h/2*（fun2（x（n）,u（n））+fun2（x（n+1）,v（k）））;

end

u（n+1）=v（k+1）;

end

plot（x,u,'o'）;

holdon

u_real=x.*（（8-3.*log（x））.^（1/3））;

plot（x,u_real,'r'）;

由图像可知：

改进的Euler法和Euler法都很接近真值。

（3）

真解

.

解：

functiondu=fun3（x,u）

du=u/（2*x）-x/（2*u^2）;

clear;

h=0.1;N=0.5/h;

x=1:

h:

1.5;

x

（1）=1;u

（1）=1;

forn=1:

N

u（n+1）=u（n）+h*fun3（x（n）,u（n））;

end

plot（x,u,'*'）;

holdon

forn=1:

N

v

（1）=u（n）+h*fun3（x（n）,u（n））;

fork=1:

6

v（k+1）=u（n）+h/2*（fun3（x（n）,u（n））+fun3（x（n+1）,v（k）））;

end

u（n+1）=v（k+1）;

end

plot（x,u,'o'）;

holdon

u_real=（4*x.^（3/2）-3*x.^2）.^（1/3）;

plot（x,u_real,'r'）;

由图像可以知道:

显然改进的Euler法要比Euler法精确度要高；

2.试用预报校正格式（1.20）解初值问题

并与Euler格式比较精度，取h=0.1。

作业要求：

写出程序，列表或用图形显示结果，并给出图或表所说明的结果。

解：

functiondu=Euler_fun2（t,u）

du=-u+t+1;

clear;

h=0.1;tend=1;N=1/h;t

（1）=0;u

（1）=1;

t=h.*（0:

N）;

forn=1:

N

u（n+1）=u（n）+h*Euler_fun2（t（n）,u（n））;

end

plot（t,u,'*'）;

holdon

forn=1:

N

u0（n+1）=u（n）+h*Euler_fun2（t（n）,u（n））;

u（n+1）=u（n）+h/2*（Euler_fun2（t（n）,u（n））+Euler_fun2（t（n+1）,u0（n+1）））;

end

plot（t,u,'o'）;

holdon

sol=dsolve（'Du=-u+t+1','u（0）=1'）;

u_real=eval（sol）;

plot（t,u_real,'r'）;

由图像可以知道:

显然预报校正格式要比Euler法精确度要高；

P37例4.1用四级四阶Runge-Kutta法计算初值问题：

u’=4tu0.5,0≤t≤2,

u（0）=1.

取h=0.1,0.5,1.精确解为u（t）=（1+t2）2

作业要求：

写出程序，列表或用图形显示结果，并给出图或表所说明的结果.

解：

functiondu=fun4（t,u）

du=4*t*u.^（1/2）;

clear;

h=0.1;N=2/h;t=0:

h:

2;t

（1）=0;u

（1）=1;

forn=1:

N

k1=fun4（t（n）,u（n））;

k2=fun4（t（n）+0.5*h,u（n）+0.5*h*k1）;

k3=fun4（t（n）+0.5*h,u（n）+0.5*h*k2）;

k4=fun4（t（n）+h,u（n）+h*k3）;

u（n+1）=u（n）+h*（k1+2*k2+2*k3+k4）/6;

end

plot（t,u,'*'）;

holdon

u_real=（1+t.^2）.^2;

plot（t,u_real,'r'）;

将上述h换为0.5后图像为：

将上述h换为1后图像为：

从上述图像来看：

第一幅图说明四级四阶Runge-Kutta法精度很高；后面两幅图说明了步长h取的越小越逼近精确值。