﻿ 140211020 邓美锋 13 - 计算数学达人 - 专，学者，数值代数，微分方程数值解

140211020 邓美锋 13

3.Use separation of variables to find solutions of the IVP given by y(0)=1 and the following differential equations:

(a)y’=t.      (d)y’=5(t^4)y.

dy/dt=t,dy=tdt,y=1/2t^2+C;

y(0)=0+C=1;

C=1;

y=1/2t^2+1.

dy/dt=5(t^4)y,(1/y)dy=(5t^4)dt;

ln|y|=t^5+C;

|y|=e^(t^5+C);

y(0)=e^C=1,C=0;

So,y=e^(t^5).

2.Plot the Euler’s Method approximate solutions for the IVPs in Exercise 3 on [0,1] for step sizes h=0.1,0.05,and 0.025,along with the exact solution.

(1)Code:

function [t,y]=euler(inter,y0,n)

t(1)=inter(1);y(1)=y0;

h=(inter(2)-inter(1))/n;

for i=1:n

t(i+1)=t(i)+h;

y(i+1)=eulerstep(t(i),y(i),h);

end

plot(t,y)

function y=eulerstep(t,y,h)

y=y+h*ydot(t,y);

function z=ydot(t,y)

z=t;

(2)Result:

>> [t,y]=euler([0 1],1,10);

>> hold on;

>> t=0:0.01:1;

>> y=(1/2)*(t.^2)+1;

>> plot(t,y)

>> [t,y]=euler([0 1],1,20);

>> hold on;

>> t=0:0.01:1;

>> y=(1/2)*(t.^2)+1;

>> plot(t,y)

>> [t,y]=euler([0 1],1,40);

>> hold on;

>> t=0:0.01:1;

>> y=(1/2)*(t.^2)+1;

>> plot(t,y)

(1)Code:

function [t,y]=euler(inter,y0,n)

t(1)=inter(1);y(1)=y0;

h=(inter(2)-inter(1))/n;

for i=1:n

t(i+1)=t(i)+h;

y(i+1)=eulerstep(t(i),y(i),h);

end

plot(t,y)

function y=eulerstep(t,y,h)

y=y+h*ydot(t,y);

function z=ydot(t,y)

z=5*(t^4)*y;

(2)Result:

>> [t,y]=euler([0 1],1,10);

>>  hold on;

>> t=0:0.01:1;

>> y=exp(t.^5);

>>  plot(t,y)

>> [t,y]=euler([0 1],1,20);

>> hold on;

>> t=0:0.01:1;

>> y=exp(t.^5);

>> plot(t,y)

>> [t,y]=euler([0 1],1,40);

>> hold on;

>> t=0:0.01:1;

>> y=exp(t.^5);

>> plot(t,y)

6.For the initial value problems in Exercise 4,make a log-log plot of the error of Euler’s Method at t=2 as a function of h=0.1×2^(-k) for 0k5.

function [t,y]=euler(inter,y0,n)

t(1)=inter(1);y(1)=y0;

h=(inter(2)-inter(1))/n;

for i=1:n

t(i+1)=t(i)+h;

y(i+1)=eulerstep(t(i),y(i),h);

end

plot(t,y)

function y=eulerstep(t,y,h)

y=y+h*ydot(t,y);

(a)(1)Code:

function z=ydot(t,y)

z=t+y;

>> [t,y]=euler([0 2],0,10);

>> t=2;error=-t-1+exp(t)-y(11)

error =

1.1973

>> [t,y]=euler([0 2],0,20);

>> t=2;error=-t-1+exp(t)-y(21)

error =

0.6616

>> [t,y]=euler([0 2],0,40);

>> t=2;error=-t-1+exp(t)-y(41)

error =

0.3491

>> [t,y]=euler([0 2],0,80);

>> t=2;error=-t-1+exp(t)-y(81)

error =

0.1795

>> [t,y]=euler([0 2],0,160);

>> t=2;error=-t-1+exp(t)-y(161)

error =

0.0910

>> [t,y]=euler([0 2],0,320);

>> t=2;error=-t-1+exp(t)-y(321)

error =

0.0458

>> h=[0.10000 0.05000 0.02500 0.01250 0.00625 0.00312];

>> error=[1.1973 0.6616 0.3491 0.1795 0.0910 0.0458];

>> loglog(h,error)

(b)(1)Code:

function z=ydot(t,y)

z=t-y;

>> [t,y]=euler([0 2],0,10);

>>  t=2;error=t-1+exp(-t)-y(11)

error =

0.0280

>> [t,y]=euler([0 2],0,20);

>> t=2;error=t-1+exp(-t)-y(21)

error =

0.0138

>> [t,y]=euler([0 2],0,40);

>> t=2;error=t-1+exp(-t)-y(41)

error =

0.0068

>> [t,y]=euler([0 2],0,80);

>> t=2;error=t-1+exp(-t)-y(81)

error =

0.0034

>> [t,y]=euler([0 2],0,160);

>> t=2;error=t-1+exp(-t)-y(161)

error =

0.0017

>> [t,y]=euler([0 2],0,320);

>> t=2;error=t-1+exp(-t)-y(321)

error =

8.4673e-004

>> h=[0.10000 0.05000 0.02500 0.01250 0.00625 0.00312];

>> error=[0.0280 0.0138 0.0068 0.0034 0.0017 8.4673e-004];

>> loglog(h,error)

(c)(1)Code:

function z=ydot(t,y)

z=4*t-2*y;

>> [t,y]=euler([0 2],0,10);

>> t=2;error=2*t-1+exp(-2*t)-y(11)

error =

0.0123

>> [t,y]=euler([0 2],0,20);

>> t=2;error=2*t-1+exp(-2*t)-y(21)

error =

0.0068

>> [t,y]=euler([0 2],0,40);

>> t=2;error=2*t-1+exp(-2*t)-y(41)

error =

0.0035

>> [t,y]=euler([0 2],0,80);

>> t=2;error=2*t-1+exp(-2*t)-y(81)

error =

0.0018

>> [t,y]=euler([0 2],0,160);

>> t=2;error=2*t-1+exp(-2*t)-y(161)

error =

9.0805e-004

>> [t,y]=euler([0 2],0,320);

>> t=2;error=2*t-1+exp(-2*t)-y(321)

error =

4.5597e-004

>> h=[0.10000 0.05000 0.02500 0.01250 0.00625 0.00312];

>> error=[0.0123 0.0068 0.0035 0.0018 9.0805e-004 4.5597e-004];

>> loglog(h,error)