﻿ 15074137杨子仪-13 chapter6 - 计算数学达人 - 专，学者，数值代数，微分方程数值解

15074137杨子仪-13 chapter6

6.1 Initial Value Problems

Exercise6.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= 5t4y

Solution

(a) The exact solution of the equation y=t can be found by using the method of separation of variables.Separate variables,and integrate,as follows:

dy=tdt

y=1/2*t^2+k

The initial condition y(0)=1,so k=1 implies y=1/2*t^2+1

(d) The exact solution of the equation y’=t can be found by using the method of separation of variables.Assuming that y0divide both sides by y,separate variables,and integrate,as follows:

dy/y=5t^4dt

ln|y|=t^5+k

|y|=e^(t^5+k)

The initial condition y(0)=1,so k=0 implies y=e^(t^5)

Computer Problems6.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.

Solution

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 z=ydot(t,y)

z=t;

function y=eulerstep(t,y,h)

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

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)

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

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 z=ydot(t,y)

z=t;

function y=eulerstep(t,y,h)

y=y+h*ydot(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)