Question

Differential Equations

1. Create a direction field for y 0 = y − y 2 .

(a) You should find any equilibrium solutions by hand and at least a few other solutions. Feel free to make a direction field with some piece of technology and share a picture of it.

(b) Find at least a few solution curves and describe the behavior of y as x → ∞, for different ranges of initial values y(x0) = y0.

(c) Use your direction field to approximate the value of y(0.5) if the initial condition is y(0) = 0.5.

(d) Use Euler’s method with h = 0.1 to approximate y(0.5) when the initial condition is y(0) = 0.5.

(e) Bonus: Find the analytical solution to this differential equation with initial condition y(0) = 0.5 and then find the exact value of y, what is your percent error from your Euler method approximation?

Answer 1

%Matlab code for direction field and plot
clear all
close all

%function for which direction field have to do
f=@(t,y) y-y.^2;
fprintf('function for which direction field have to make f(y)=')
disp(f)
%Plotting Direction field
figure(1)
dirrfield(f,0:.25:5,-2:.25:4)
hold on
t_ini=0;
t_end=5;
%plot for different initial velocity
for v0=-2:0.5:4
[ts,ys] = ode45(f,[t_ini,t_end],v0); plot(ts,ys,'b','Linewidth',2)
end

hold off
title('phase potrait plot')
xlabel('time in sec.')
ylabel('y(t) ')

syms y(t) y0
eqn = diff(y,t) == y-y.^2;
cond = [y(0)==y0];
ySol(t) = dsolve(eqn,cond);
fprintf('Exact solution for initial condition y(0)=y0, is ')
disp(ySol)

fprintf('For y(0)=0.5, y(0.5)=0.62 \n')
fprintf('Solution using Euler method\n')
[y1_result,t_result] = euler_method(f,0.5,0,0.5,0.1);
%Exact solution for y(0)=0.5
syms y(t)
eqn = diff(y,t) == y-y.^2;
cond = [y(0)==0.5];
yS(t) = dsolve(eqn,cond);
fprintf('Exact solution for initial condition y(0)=0.5, is ')
disp(yS)

for i=1:6
    err=abs(y1_result(i)-yS(t_result(i)))/yS(t_result(i));
    fprintf('\tExact solution at t=%2.2f is %f with error=%f\n',t_result(i),y1_result(i),err)
end

%%Matlab function for direction filed
function dirrfield(f,tval,yval)
% dirfield(f, t1:dt:t2, y1:dy:y2)
%
%   plot direction field for first order ODE y' = f(t,y)
%   using t-values from t1 to t2 with spacing of dt
%   using y-values from y1 to t2 with spacing of dy
%
%   f is an @ function, or an inline function,
%     or the name of an m-file with quotes.
%
% Example: y' = -y^2 + t
%   Show direction field for t in [-1,3], y in [-2,2], use
%   spacing of .2 for both t and y:
%
%   f = @(t,y) -y^2+t
%   dirfield(f, -1:.2:3, -2:.2:2)

[tm,ym]=meshgrid(tval,yval);
dt = tval(2) - tval(1);
dy = yval(2) - yval(1);
fv = vectorize(f);
if isa(f,'function_handle')
fv = eval(fv);
end
yp=feval(fv,tm,ym);
s = 1./max(1/dt,abs(yp)./dy)*0.75;
h = ishold;
quiver(tval,yval,s,s.*yp,0,'r'); hold on;
%quiver(tval,yval,-s,-s.*yp,0,'r');
if h
hold on
else
hold off
end
axis([tval(1)-dt/2,tval(end)+dt/2,yval(1)-dy/2,yval(end)+dy/2])
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%Matlab code for Euler's forward
function [y1_result,t_result] = euler_method(f,y0,t0,tend,h)
%function for Euler equation solution
  
    %all step size
    N=(tend-t0)/h;
    %Initial values

    %t end values
    tn=t0:h:tend;
    % Euler steps
    y1_result(1)=y0;
    t_result(1)=t0;
    fprintf('\nFor step size %2.2f\n',h)
    for i=1:length(tn)-1
        t_result(i+1)= t_result(i)+h;
        y1_result(i+1)=y1_result(i)+h*double(f(t_result(i),y1_result(i)));
        fprintf('\tAt t=%2.2f y(%2.2f)=%f\n',t_result(i+1),t_result(i+1),y1_result(i+1))
    end
end

%%%%%%%%%%%%%% End of Code %%%%%%%%%%%%%%