﻿ 140231009+杨蕾+13+Initial Value Problems - 计算数学达人 - 专，学者，数值代数，微分方程数值解

### 140231009+杨蕾+13+Initial Value Problems

P291  Exercises 6.1

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;

(a)

dy/dt=t,  dy=t*dt,  ∫dy=∫t*dt,

y=1/2*t^2+c,  1=y(0)=c,

Therefore,y=1/2*t^2+1.

(d)

dy/dt=5*t^4*y,  dy/y=5*t^4*dt,  ∫dy/y=∫5*t^4*dt,

ln|y|=t*5+c1,  |y|=e^(t*5+c1),  1=y(0)=e^c1,c1=0,

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

P292  Computer Problems 6.1

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.

(a)y'=t;

function [t,y]=Eu(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+0*y;

Results:

(1)h=0.1:

>> x=0:0.1:1;

>> y=1/2*x.^2+1;

>> plot(x,y)

>> hold on

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

t =

0    0.1000    0.2000    0.3000    0.4000    0.5000    0.6000    0.7000    0.8000    0.9000    1.0000

y =

1.0000    1.0000    1.0100    1.0300    1.0600    1.1000    1.1500    1.2100    1.2800    1.3600    1.4500

Picture:

(2)h=0.05:

>> x=0:0.05:1;

>> y=1/2*x.^2+1;

>> plot(x,y)

>> hold on

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

t =

1 14

0    0.0500    0.1000    0.1500    0.2000    0.2500    0.3000    0.3500    0.4000    0.4500    0.5000    0.5500    0.6000    0.6500

15 21

0.7000    0.7500    0.8000    0.8500    0.9000    0.9500    1.0000

y =

1 14

1.0000    1.0000    1.0025    1.0075    1.0150    1.0250    1.0375    1.0525    1.0700    1.0900    1.1125    1.1375    1.1650    1.1950

15 21

1.2275    1.2625    1.3000    1.3400    1.3825    1.4275    1.4750

Picture:

(3)h=0.025:

>> x=0:0.025:1;

>> y=1/2*x.^2+1;

>> plot(x,y)

>> hold on

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

t =

1 14

0    0.0250    0.0500    0.0750    0.1000    0.1250    0.1500    0.1750    0.2000    0.2250    0.2500    0.2750    0.3000    0.3250

15 28

0.3500    0.3750    0.4000    0.4250    0.4500    0.4750    0.5000    0.5250    0.5500    0.5750    0.6000    0.6250    0.6500    0.6750

29 41

0.7000    0.7250    0.7500    0.7750    0.8000    0.8250    0.8500    0.8750    0.9000    0.9250    0.9500    0.9750    1.0000

y =

1 14

1.0000    1.0000    1.0006    1.0019    1.0038    1.0063    1.0094    1.0131    1.0175    1.0225    1.0281    1.0344    1.0413    1.0488

15 28

1.0569    1.0656    1.0750    1.0850    1.0956    1.1069    1.1188    1.1313    1.1444    1.1581    1.1725    1.1875    1.2031    1.2194

29 41

1.2363    1.2538    1.2719    1.2906    1.3100    1.3300    1.3506    1.3719    1.3938    1.4163    1.4394    1.4631    1.4875

Picture:

(d)y'=5*t^4*y;

function [t,y]=Eu(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;

Results:

(1)h=0.1:

>> x=0:0.1:1;

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

>> plot(x,y)

>> hold on

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

t =

0    0.1000    0.2000    0.3000    0.4000    0.5000    0.6000    0.7000    0.8000    0.9000    1.0000

y =

1.0000    1.0000    1.0001    1.0009    1.0049    1.0178    1.0496    1.1176    1.2517    1.5081    2.0028

Pictures:

(2)h=0.05:

>> x=0:0.05:1;

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

>> plot(x,y)

>> hold on

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

t =

1 10

0    0.0500    0.1000    0.1500    0.2000    0.2500    0.3000    0.3500    0.4000    0.4500

11 20

0.5000    0.5500    0.6000    0.6500    0.7000    0.7500    0.8000    0.8500    0.9000    0.9500

21

1.0000

y =

1 10

1.0000    1.0000    1.0000    1.0000    1.0002    1.0006    1.0015    1.0036    1.0073    1.0138

11 20

1.0242    1.0402    1.0640    1.0984    1.1475    1.2163    1.3125    1.4469    1.6358    1.9041

21

2.2918

Picture:

(3)h=0.025:

>> x=0:0.025:1;

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

>> plot(x,y)

>> hold on

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

t =

1 10

0    0.0250    0.0500    0.0750    0.1000    0.1250    0.1500    0.1750    0.2000    0.2250

11 20

0.2500    0.2750    0.3000    0.3250    0.3500    0.3750    0.4000    0.4250    0.4500    0.4750

21 30

0.5000    0.5250    0.5500    0.5750    0.6000    0.6250    0.6500    0.6750    0.7000    0.7250

31 40

0.7500    0.7750    0.8000    0.8250    0.8500    0.8750    0.9000    0.9250    0.9500    0.9750

41

1.0000

y =

1 10

1.0000    1.0000    1.0000    1.0000    1.0000    1.0000    1.0000    1.0001    1.0002    1.0004

11 20

1.0007    1.0012    1.0020    1.0030    1.0044    1.0063    1.0087    1.0120    1.0161    1.0213

21 30

1.0278    1.0358    1.0457    1.0576    1.0721    1.0894    1.1102    1.1350    1.1645    1.1994

31 40

1.2408    1.2899    1.3481    1.4171    1.4991    1.5970    1.7140    1.8545    2.0243    2.2303

41

2.4823

Picture:

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 0≤k≤5.

(a)y'=t+y;

function Plot(inter,y0)

t(1)=inter(1);

y(1)=y0;

for k=0:1:5

h(k+1)=0.1*2.^(-k);

n=2/h(k+1);

for i=1:n

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

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

end

e(k+1)=y(n+1)-(exp(2)-3);

end

loglog(h,e,'.')

function y=eulerstep(t,y,h)

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

function z=ydot(t,y)

z=t+y;

Results:

>> Plot([0 2],0)

Picture:

(b)y'=t-y;

Code:

function Plot(inter,y0)

t(1)=inter(1);

y(1)=y0;

for k=0:1:5

h(k+1)=0.1*2.^(-k);

n=2/h(k+1);

for i=1:n

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

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

end

e(k+1)=y(n+1)-(exp(-2)+1);

end

loglog(h,e,'.')

function y=eulerstep(t,y,h)

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

function z=ydot(t,y)

z=t-y;

Results:

>> Plot([0 2],0)

Picture:

(c)y'=4t-2y;

Code:

function Plot(inter,y0)

t(1)=inter(1);

y(1)=y0;

for k=0:1:5

h(k+1)=0.1*2.^(-k);

n=2/h(k+1);

for i=1:n

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

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

end

e(k+1)=y(n+1)-(exp(-4)+3);

end

loglog(h,e,'.')

function y=eulerstep(t,y,h)

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

function z=ydot(t,y)

z=4*t-2*y;

Results:

>> Plot([0 2],0)

Pictures: