﻿ 15074143+欧阳章敏+13+6.1 - 计算数学达人 - 专，学者，数值代数，微分方程数值解

### 15074143+欧阳章敏+13+6.1

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)

dy/dt=t

dy=t dt

y=1/2*t^2+c

According to y(0)=1 ,hence, c=1.So, y=1/2*t^2+1.

(d)

dy/dt=5*t^4*y

1/y dy=5*t^4 dt

log y= t^5 + c

y=c*e^t^5

According to y(0)=1 ,hence, c=1.So, 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 x=euler(f,x0,y0,xn,h)

n=(xn-x0)/h;

x=zeros(1,n+1);

y=zeros(1,n+1);

x(1)=x0;

y(1)=x0;

for n=1:n

x(n+1)=x(n)+h;

y(n+1)=y(n)+h*feval(f,x(n),y(n));

end

plot(x,y)

function z=li9_1fun(x,y)

z=x;

x0=0;y0=1;xn=1;h=0.1;
x=euler('li9_1fun',x0,y0,xn,h);
>> hold on
>> x=0:0.1:1;
>> y=dsolve('Dy=x','y(0)=1','x');
>> ezplot(y)

x0=0;y0=1;xn=1;h=0.05;
>> x=euler('li9_1fun',x0,y0,xn,h);
>>  hold on
>> x=0:0.05:1;
>> y=dsolve('Dy=x','y(0)=1','x');
>> ezplot(y)

x0=0;y0=1;xn=1;h=0.025;
x=euler('li9_1fun',x0,y0,xn,h);
hold on
x=0:0.025:1;
y=dsolve('Dy=x','y(0)=1','x');
ezplot(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.Answer:

(a):

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=1/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(3)-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)

Pictures:

(b):

Code:

function Plot2(inter,y0)
t(1)=inter(1);
y(1)=y0;
for k=0:1:5
h(k+1)=0.1*2.^(-k);
n=1/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:

>> Plot2([0 2],0)

Pictures:

(c):
Code:

function Plot3(inter,y0)
t(1)=inter(1);
y(1)=y0;
for k=0:1:5
h(k+1)=0.1*2.^(-k);
n=1/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:

>> Plot3([0 2],0)

Pictures: