﻿ 15074139胡雨薇-13-IVP - 计算数学达人 - 专，学者，数值代数，微分方程数值解

### 15074139胡雨薇-13-IVP

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：
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 y≠0，divide 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:
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;

When h=0.1

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

>> hold on;

>> t=0:0.01:1;

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

>> plot(t,y)

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

>> hold on;

>> t=0:0.01:1;

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

>> plot(t,y)

When h=0.025

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

>> hold on;

>> t=0:0.01:1;

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

>> 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.
Solution:
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);

(a)(1)

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)
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)