现在时间是:
当前位置:首 页 >> 数值分析(英)>> 教学区>> 文章列表

140231009+杨蕾+13+Initial Value Problems

作者:   发布时间:2017-06-14 13:07:37   浏览次数:31

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;

Answer:

(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;

Answer & Code:

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;

Answer & Code:

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;

Answer & 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)-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:

 







Copyright ©2017    计算数学达人 All Right Reserved.

Powered by LingCms 版权所有 © 2005-2017 Lingd.Net.粤ICP备16125321号 -4