Differential Equations: Homework on Phase Curves
Improved Euler Method
(T, Y, Z)
1.0e+006 *
0.0000 0.0000 -0.0000
0.0000 0.0000 -0.0000
0.0000 0.0000 -0.0000
0.0000 -0.0000 -0.0000
0.0000 -0.0000 -0.0000
0.0000 -0.0000 -0.0000
0.0000 -0.0000 -0.0000
0.0000 -0.0001 -0.0003
0.0000 -0.0083 -0.0414
0.0000 -1.2292 -6.1462
Actual Values =
(T, Y)
0.9990 0.0068
1.9990 0.0000
2.9990 0.0000
3.9990 0.0000
4.9990 0.0000
5.9990 0.0000
6.9990 0.0000
7.9990 0.0000
8.9990 0.0000
9.9990 0.0000
This graph shows Z on the vertical axis and Y on the horizontal. From a comparison
of the actual values of Y vs. the Improved Euler approximation, we can see that
the error is quite large. The general solution of the equation is of the form
a*exp(-5t) + b*exp(5t). For our initial conditions b is zero, but just as the
solution curve deviates slightly from the correct curve for our given initial
conditions, we can see that the b*exp(5t) part of the solution starts to contribute,
and error builds up as the approximation proceeds.
Code:
%
Title: Phase Curves with Improved Euler
%
Author: Bobby Rohrkemper
%
Kalamazoo College
%
Date Created: February 17, 2002
%
draws a direction field
%
(function, default window
dirfldyz('(25*y)/z')
% clears all x, y, p, and actual values
t = []; y = []; z = []; t1 = []; y1 = []; z1 = []; actualY = []; actualY1
= [];
%
x1 and y1 are set to their initial conditions
t(1) = 0; y(1) = 1; z(1) = -5; actualY(1) = 1;
%
sets the timestep and the number of steps
h = .001; numSteps = 10000;
%
loops through a number of timesteps and
%
finds a solution using the improved Euler method
k = 1; % used in if statement
for
i = 1 : numSteps, % number of timesteps
t(i+1) = t(i) + h; % computes new t
fOld = 25 * y(i); % computes the old value
%y(i+1) = y(i) + h * z(i);
%z(i+1) = z(i) + h * fOld;
pY = y(i) + h * z(i);
fNew = 25 * pY ;
pZ = z(i) + h * fOld;
y(i+1) = y(i) + (h/2) * ( z(i) + pZ );
z(i+1) = z(i) + (h/2) * ( fOld + fNew );
% compute the actual y values
actualY(i+1) = exp( -5 * t(i+1) );
% stores every nth value computed for t, y, z
% (in this case, n is 1000)
if( mod(i,1000) == 0 )
t1(k) = t(i);
y1(k) = y(i);
z1(k) = z(i);
actualY1(k) = actualY(i);
k = k + 1;
;end
% compute the actual y values
actualY(i+1) = exp( -5 * t(i+1) );
;end
% prints out a table of x values
[ t1' y1' z1' ]
[ t1' actualY1' ]
% plots the computed solution in a separate window with its direction field
plot(y,
z, 'g')
19.
(a.)
Z’ = - y3
+ y
(b.)
Z’ = -y3
+ 2y
(c.) Z’ = -2y3 + y
Code:
Here is the code for part
c:
% Author: Bobby Rohrkemper
% Kalamazoo
College
% Date
Created: February 17, 2002
% draws
a direction field
% (function,
default window
dirfldyz('(-y^3 + 2*y)/z')
% clears all x, y, p, and actual values
t = []; y = []; z = []; t1 = []; y1 = []; z1 = []; actualY = []; actualY1
= [];
% x1 and
y1 are set to their initial conditions
t(1) = 0; y(1) = 1.5; z(1) = 0; actualY(1) = 0;
for j
= 0 : 10;
y(1) = j/2;
% sets the timestep and the number of steps
h = .001; numSteps = 40000;
% loops through a number of timesteps and
% finds a solution using the improved Euler method
k = 1; % used in if statement
for i = 1 : numSteps, % number of timesteps
t(i+1) = t(i) + h; % computes new t
fOld = ( -y(i)^3 + 2*y(i) ); % computes the old value
%y(i+1) = y(i) + h * z(i);
%z(i+1) = z(i) + h * fOld;
pY = y(i) + h * z(i);
fNew = (-pY^3 + 2*pY) ;
pZ = z(i) + h * fOld;
y(i+1) = y(i) + (h/2) * ( z(i) + pZ );
z(i+1) = z(i) + (h/2) * ( fOld + fNew );
;end
% plots
the computed solution in a separate window with its direction field
plot(y, z, 'g')
;end
Here is a phase graph of the equation for a hard spring oscillator. We can
see that there is a portion of the graph to the right of the Z axis in which
there is a constant amplitude of (3/10).
Looking at a close up of
this section to the right of the Z axis in a window with vertical limits of
-.3 and .3, we can see that the amplitude is almost exactly .3.
We can realize that the amplitude of this portion comes from the (3/10)cos(t)
term based on what we know about oscillatory system equations. If we re-arrange
the equation so that the cos(t) term is the only one on the right side of the
equation, we can see this term as the driving force. We know that the system
will eventually move in a way completely defined by the driving force after
enough time has passed. So, we have identified a phase curve in the phase portrait
that represents a periodic path. The switch to periodic behavior seems to happen
after y becomes positive, and it stays that way for an indefinite amount of
time.
Source: 2ORDPLOT java applet
code:
% Title:
Phase Curves with Improved Euler
% Author:
Bobby Rohrkemper
% Kalamazoo
College
% Date
Created: February 17, 2002
% draws
a direction field
% (function,
default window
dirfldyz('(exp(-2*y)-exp(-y))/z')
% clears all x, y, p, and actual values
t = []; y = []; z = []; t1 = []; y1 = []; z1 = []; actualY = []; actualY1
= [];
% x1 and
y1 are set to their initial conditions
t(1) = 0; y(1) = 1.5; z(1) = 0; actualY(1) = 0;
for j
= 0 : 10;
y(1) = j/2;
% sets the timestep and the number of steps
h = .001; numSteps = 30000;
% loops through a number of timesteps and
% finds a solution using the improved Euler method
k = 1; % used in if statement
for i = 1 : numSteps, % number of timesteps
t(i+1) = t(i) + h; % computes new t
fOld = ( (exp(-2*y(i))-exp(-y(i))) ); % computes the old
value
%y(i+1) = y(i) + h * z(i);
%z(i+1) = z(i) + h * fOld;
pY = y(i) + h * z(i);
fNew = (exp(-2*pY)-exp(-pY)) ;
pZ = z(i) + h * fOld;
y(i+1) = y(i) + (h/2) * ( z(i) + pZ );
z(i+1) = z(i) + (h/2) * ( fOld + fNew );
;end
% plots
the computed solution in a separate window with its direction field
plot(y, z, 'g')
;end
|
