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

140211005施莎莎—作业13

作者:施莎莎   发布时间:2017-06-19 08:19:55   浏览次数:39

6.1 Exercises

 

Problems 3.
Use separation of variables to find solutions of the IVP given by y(0)=1 and the 
following differential equations:
 (a)y'=t;
Answer:
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.
 
 (b)y'=5*t^4*y;
Answer:
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,
y=e^(t^5).
 
6.1 Computer Problems.
Problem 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;
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:
 
 (b)y'=5*t^4*y;
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:
 
Problem 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;
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