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

15074139胡雨薇-13-IVP

作者:15074139   发布时间:2017-06-14 13:19:56   浏览次数:42

 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)
 






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

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