微分方程作业

P10习题

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

解：function du=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); for n=1:N

u(n+1)=u(n)+h*Euler_fun1(t(n),u(n)); end

plot(t,u,'*');hold on for n=1:N

v(1)=u(n)+h*Euler_fun1(t(n),u(n)); for k=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；并比较两个算法的精度。

解：

function du=fun31(y) du=y;

function dy=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); for n=1:N

u(n+1)=u(n)+h*y(n); y(n+1)=y(n)+h*(-u(n)); end

plot(t,u,'*');hold on for n=1:N

v(1)=u(n)+h*fun31(y(n)); w(1)=y(n)+h*fun32(u(n));

for k=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.取步长 h?0.1 ，分别用Euler 法和改进的Euler 法求下列初值问题的解，并与真解

相比较.

2x??u'?u?,0?x?1,（1）? u??u(0)?1,真解 u(x)?1?2x； 解：

function du=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法——%

for n=1:N

u(n+1)=u(n)+h*fun1(x(n),u(n)); end

plot(x,u,'*'); hold on

%——改进的Eluer法——%

for n=1:N

v(1)=u(n)+h*fun1(x(n),u(n)); for k=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'); hold on

%——真解——%

u_real=sqrt(1+2*x); plot(x,u_real,'r');

由图像可以知道:

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

?ux2?u'??2,1?x?2,（2）? xu??u(1)?2,真解 u(x)?x(8?3lnx)； 解：

function du=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; for n=1:N

u(n+1)=u(n)+h*fun2(x(n),u(n)); end

plot(x,u,'*'); hold on for n=1:N

v(1)=u(n)+h*fun2(x(n),u(n)); for k=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'); hold on

u_real=x.*((8-3.*log(x)).^(1/3)); plot(x,u_real,'r');

13

由图像可知：

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

ux???2,1?x?1.5,（3）?u'?2x2u

?u(1)?1,?真解 u(x)?(4x?3x). 解：function du=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; for n=1:N

u(n+1)=u(n)+h*fun3(x(n),u(n)); end

plot(x,u,'*'); hold on for n=1:N

v(1)=u(n)+h*fun3(x(n),u(n)); for k=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'); hold on

u_real=(4*x.^(3/2)-3*x.^2).^(1/3); plot(x,u_real,'r');

32123

由图像可以知道:

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

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

t??0,1??u???u?t?1??u|t?0?1并与Euler格式比较精度，取h=0.1。作业要求：写出程序，列表或用图形显示结果，并给出图或表所说明的结果。 解：

function du=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); for n=1:N

u(n+1)=u(n)+h*Euler_fun2(t(n),u(n)); end

plot(t,u,'*'); hold on for n=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'); hold on

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

作业要求：写出程序，列表或用图形显示结果，并给出图或表所说明的结果. 解：

function du=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; for n=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,'*');

hold on

u_real=(1+t.^2).^2; plot(t,u_real,'r');

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

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

从上述图像来看：

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

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