Final Exam - Answers

(Nearly all of these were taken directly from the textbook, 
	Boyce-DiPrima 9th "Elementary Differential Equations",
	as noted below by each problem. )

Differential Equations - Spring 2009 - Marlboro College

Please do read the instructions first.  Good luck.  - Jim M

instructions

Your mission is to use the six differential equations below 
to showcase your understanding of the material we've covered this term, namely
	* exact solutions, including separation of variables, integrating factors, and so on,
	* changing the physical dimensions to find a convenient form of the equation,
	* graphical approaches, including direction fields and solution curves, and phase space
	* numerical methods such as runge-kutta and heun,
	* initial value problems and integration constants
	* system of first order equations,
	* the use of eigenvalues and eigenvectors,
	* critical points and the qualitative long-term behavior of the solutions,
	* stability and/or chaos in particular circumstances.

All greek letters (e.g. α,β,γ,Λ,...) are constants, 
while roman letters (e.g. x, y, z, u, v, t, ...) are variables.

While you shouldn't try to include all of these techniques in every equation,
most of them will be appropriate for at least one of the equations.  And many
of the equations will naturally use several of these ideas.  I'm expecting exact solutions 
for some, graphs of diverse sorts for several, numerical solutions from a 
particular initial conditions for others, and so on.  You're welcome to use more 
than one approach on the same problem, for example adding a numeric solution to one 
that you can do analytically, as a way to look at the accuracy of the numerical method.

In each case I've simply given an equation, and not specified a particular initial condition;
however, I encourage you to choose an initial condition (or several) and work from there.
You may also want to choose specific values for any constants - or even turn numbers
like 4 into constants like β and discuss what happens for similar equations.  Whatever
you do, be clear about your choices and assumptions.

You have one day to do the exam, which is due by email, posted on the assigments 
page, and/or tucked under my office door by midnight Friday May 8.  You're welcome 
to use any format you like: Mathematica, Excel, python, by hand, or whatever.

Sources such as our text or the 'net are OK; however, if so you must cite them.
Remember that you're being evaluated on your understanding of the material,
as I can see it in what you write - not just "the right answer" - so be clear in 
your presentation.  

Asking other people for help is not allowed.

If you have questions, you can email me or phone at either 247-0857 (home) or 258-9255 (office).

code

Rather than intersperse algorithms with the answers,
I've put here the numerical routines that I use later.

To keep this simple, I'm only using two : 
one to plot direction fields, and one to find solutions from an initial (x0, y0) point.

Both these routines are from my nonlinear1.nb notebook, from my April 23 notes.
Heun's algorithm could be used instead of RK; the result would be similar but with either 
less accuracy or a smaller time step.

In each case, the equation(s) must be put in the form of two first order equations
(i.e. dx/dt=..., dy/dt=...) to use rungekutta2[], and a slope dy/dx must be found to use 
the directionField[].  

Problems 4, 5, and 6 below can all be put into a form that suits these routines.

(* Return a plot of a direction field for dydx[x,y]. *)
directionField[dydx_,xmin_,xmax_,ymin_,ymax_,xLabel_,yLabel_,title_]:=
VectorPlot[
{1.0/Sqrt[1+((dydx[x,y])^2)],1.0*dydx[x,y]/Sqrt[1+((dydx[x,y])^2)]},
{x,xmin,xmax},{y,ymin,ymax},
VectorStyle->"Segment",VectorScale->{0.02,Automatic,Automatic},
FrameLabel->{xLabel,yLabel},
PlotLabel->title
];

(* Return a list {{x1,y1},{x2,y2},...} *)(* of points from the solution of    *)
(* the equations dx/dt = f(x,y), dy/dt=g(x,y)   *)
(* with a 2D runge kutta 4th order,     *) 
(* starting at {t0,x0,y0}, with dt=h as the step size. *)
rungekutta2[fx_,fy_,{t0_,x0_,y0_},tEnd_, h_]:=Block[{t,x,y,k1x,k2x,k3x,k4x,k1y,k2y,k3y,k4y,points},
t=1.0*t0;
y=1.0*y0;
x=1.0*x0;
points={{x,y}};
While[t<tEnd,
k1x=fx[x,y];
k1y=fy[x,y];
k2x=fx[x+k1x h/2,y+k1y h/2];
k2y=fy[x+k1x h/2,y+k1y h/2];
k3x=fx[x+k2x h/2,y+k2y h/2];
k3y=fy[x+k2x h/2,y+k2y h/2];
k4x=fx[x+k3x h,y+k3x h];
k4y=fy[x+k3y h,y+k3y h];
t = t+h;
x=x+(k1x+2 k2x+2k3x+k4x) h/6;
y=y+(k1y+2 k2y+2k3y+k4y) h/6;
points = Append[points,{x,y}];
];
Return[points];
];

1.

((dy)/(dx))+(e^(2x))+y-1=0

1. answer

This can be done exactly with the "integrating factor" trick;
it's similar to the problems on page 39 of our text.

In the form

	y+((dy)/(dx))=1-(e^(2x))
we want the left to be a total derivative of something.  

Guessing from the "y" and "dy/dx" that there's a product of y and something, we see that
	(d/(dx))(y f(x)) = y ((df)/(dx))+((dy)/(dx))f
which works if 
	f(x)=((df)/(dx))
which is true for
	f(x)=(e^x)
So we have
	y+((dy)/(dx))=1-(e^(2x))
	((ye)^x)+((dy)/(dx))(e^x)=(e^x)(1-(e^(2x)))
	(d/(dx))(((ye)^x))=(e^x)-(e^( 3x))
And then move (dx) over to the right and integrate both sides to get
	((ye)^x)=(e^x)-(e^(3x))/3+C
or
	y=1-(1/3)(e^(2x))+C (e^(-x))		EXACT ANALYTIC ANSWER
where C is an arbitrary constant.

Check this answer by having Mathematica do the analytic differentiation :

y=1-Exp[2x]/3+C Exp[-x];
dydx=D[y,x];
dydx+Exp[2x]+y-1==0

Other things that could be done here :
	* solve for C given an initial ((Subscript[x, 0]),(Subscript[y, 0]))
	* plot some solutions for various C's
	* solve numerically and compare with exact answer,
	   using a 1D Heun or RK routine.

2.

((dy)/(dx))+α (x^2)y=0

2. answer

This one can be solved quickly by separation of variables an integration.
It's like the questions on page 49 of the text.

However, it's also the best example to use to show how it can be possible
to get rid of constants by changing the scale of the problem, in other words 
by deciding what physical units to use.

Let
	y = A u
x = B v
where A and B are constants, and (u,v) are the new variables.  
Then the equation becomes
	((d(A u))/(d(B v)))+α (((B v))^2)(A u)=0
or
	(A/B)(((du)/(dv))+α (B^3) (v^2)u)=0

We see that changing y doesn't help any; the constant A just floats out to the left,
because the original equation is linear in y.

However, we can choose B to eliminate α :
	B = (α^(1/3))
Then the equation to solve is
	((du)/(dv))+ (v^2)u=0
and we see that we don't need to consider solutions with different values of α; 
they're all the same except for a scaling of x.

Separating and integrating gives
	((du)/u)=-(v^2)dv
or
	ln(u)=-((v^3)/3)+K
or
	u=C (e^(-(v^3)/3))			EXACT ANALYTIC ANSWER
Check :

u=C Exp[-(v^3)/3];
dudv=D[u,v];
dudv + (v^2)u==0

Like problem 1, it could also be OK to show numeric solutions to this one.

3.

(((d^2)y)/(d(x^2)))-cos(x)-((dy)/(dx))-y=0

3. answer

This was supposed to be an example of the forces, damped, harmonic oscillator,
which we did in some detail.  

But I put the wrong sign in on the spring force term, (-y) rather than (y),
which means to do it right you'd need to adapt the approaches I took in class.

The notes that are similar are from the Feb 26 "oscillator2.nb" notes.

The big ideas are :
	* Guess a form of the long-term (Subscript[y, long])(x) solution 
	   in terms of exponential and trig functions
	* Solve for the constants in that guess by plugging in.
	* Find two more transient solutions (Subscript[y, 1])(x), (Subscript[y, 2])(x)  
	   to the homogeneous equation, without cos(x).
	* The general solution is then
	   	A (Subscript[y, 1])+B (Subscript[y, 2])+(Subscript[y, long])
	   where A and B are chosen to match the initial conditions.

The long term solution can be written in the form
	y = A cos(x) + B sin(x)
Plugging in and solving gives 
	A=-2/5 
B=-1/5
Check :

yy=-(2/5)Cos[x]+-(1/5)Sin[x];
dyy=D[y,x];
d2yy=D[dyy,x];
d2yy-Cos[x]-dyy-y==0

The homogeneous solutions both solve the equation 
	(((d^2)y)/(d(x^2)))-((dy)/(dx))-y=0
which is an exponential growth/decay situation.  Guessing
	y=C ((e )^(k x))
gives a solution if 
	(k^2)-k-1==0
which works for 
	(Subscript[k, ±])=(((1±(Sqrt[5])))/2)

So the overall general solution is
	y = -(2/5)cos(x)-(1/5)sin(x)+(Subscript[C, 1])(e^((Subscript[k, 1])x))+(Subscript[C, 2])(e^((Subscript[k, 2])x))		EXACT  GENERAL SOL'N

Because one of the two homogeneous solutions has a positive exponential exponent, except in special initial conditions, that one will dominate in the long run.  

This one was harder than I expected because of a typo on my part ... sorry about that.

Since this one was a 2nd order equation, another way to proceed would have been 
to convert to a system of two 1st order equations, and use a numerical approach 
like I do in the next few problems.  The initial conditions must specifiy two values,
typically y and dy/dx at x=0.

4.

(((d^2)y)/(d(x^2)))-(1-(y^2))((dy)/(dx))+y=0

4. answer

This is the "van der Pol" equation, which is discussed in some detail 
on page 559 of the text.  The solutions for this version are plotted 
on the figures on page 561.

The simplest approach is to turn it into a system of two equations, 
and in that form look it numerically.  All the solutions head towards a cyclic
limit cycle.  My nonlinear1.nb notes from Apr 23 would be a good
choice of numercial algorithms, like what I do in the text problem.

An analytic approach is to linearize it at the critical point (0,0)
and look at the eigens.

Here's the numerical approach, writing these as a system of first order equations,
defining z=dy/dx.  This gives
	((dy)/(dx))=z
((dz)/(dx))=-y+(1-(y^2))z

Using the same sorts of numerical tools we've been working with for some time now :

dydx4[y_,z_]:=z;
dzdx4[y_,z_]:=-y+(1-(y^2))z;
slope4[y_,z_]:=(-y+(1-(y^2))z)/z;

dt=0.01;
t0=0.0;
tEnd=20.0;
curv1=rungekutta2[dydx4,dzdx4,{t0,0.5,0.2},tEnd,dt];
curv2=rungekutta2[dydx4,dzdx4,{t0,-3,2},tEnd,dt];
curv3=rungekutta2[dydx4,dzdx4,{t0,3,-4},tEnd,dt];

Show[
directionField[slope4, -5,5, -5,5,"y","z","problem 4"],
ListPlot[{curv1,curv2,curv3},
PlotStyle->Red,Joined->True]
]

Cuves starting near (0,0) spiral outward and land on the cyclic heavy red line
in the figure above; those further spiral in (or "crash" inwards towards z=0),
ending up on the same curve.  

This is an example of a differential equation with a "limit cycle".  
In fact, this one only has the one limit cycle; all solutions converge 
on the same periodic behavior, regardless of where they start.

Ignoring the z=dy/dt variable and looking at just y(t), it looks like this :
(curv1 is a list of {{y1,z1}, {y2,z2}, ...}, so curve1[[i]][[1]] is the i'th y value.)

yt=Table[
{dt*i,curv1[[i]][[1]]},
{i,Range[Length[curv1]]}
];

ListPlot[yt,Joined->True,AxesLabel->{"t","y"}]

5.

((dx)/(dt))+x+y+1=0
((dy)/(dt))-2x+y-5=0

5. answer

This is problem number 15 on page 495; 
A synopsis of the answer is given on page 617, near the bottom.

An exact linear system of two equations, here I wanted you to find the eigens 
of the matrix and talk about stability.  

As a vector equation in the form we're used to dealing with, this would be

	(d/(dt))(({{x}, {y}}))=(({{-1, -1}, {2, -1}}))(({{x}, {y}}))+(({{-1}, {5}}))

which has its critical point where dx/dt=dydt=0, at (x=2,y=-1).
To see this, the matrix and its eigens are

matrix={{-1,-1},{2,-1}};
Eigenvalues[matrix]

The derivatives are zero (which can be done by matrix methods, or by setting the two equations to zero and solving for x and y) at

Inverse[matrix] . {-1,5}

Of course, finding these by hand - or explaining how - would be worth some brownie points

Because these are complex with negative real parts, the solutions spiral inwards 
around the critical point.

The general analytic solution can be written using the eigenvectors and exponentials of these
two eigenvalues (Subscript[λ, 1]), (Subscript[λ, 2]) :

	(({{y}, {x}}))=(Subscript[C, 1])(({{vec1_x}, {vec1_y}}))(e^((Subscript[λ, 1])t))+(Subscript[C, 2])(({{vec2_x}, {vec2_y}}))(e^((Subscript[λ, 2])t))

Then given some particular initial value, one can solve for specific values for the C's.

Using the directionfield and runge-kutta numerical approach, 
the direction field and some solutions can be plotted, showing 
the inward spiral behavior and the critical point at (x,y)=(-2,1)

dxdt[x_,y_]:=-x-y-1;
dydt[x_,y_]:=2x-y+5;
slope[x_,y_]:=(2x-y+5)/(-x-y-1);

dt=0.01;
t0=0.0;
tEnd=10.0;
curve1=rungekutta2[dxdt,dydt,{t0,8,0},tEnd,dt];
curve2=rungekutta2[dxdt,dydt,{t0,5,-5},tEnd,dt];
curve3=rungekutta2[dxdt,dydt,{t0,-8,0},tEnd,dt];

Show[
directionField[slope, -10,10, -10,10,"x","y","problem 5"],
ListPlot[{curve1,curve2,curve3},
PlotStyle->Red,Joined->True]
]

6.

((dx)/(dt))-(1-(((x+y))/2)) x=0
((dy)/(dt))-(2x-1) (y/4)=0

6. answer

This is problem 3 on page 540 of the text.
A brief synopsis of the "answer" (without explanation)
is at the top of page 622.

The context can be seen in the section just before that:
it's a population model, with several critical points.

The conceptually simplest approach is to look at it numerically; plot direction fields and 
some curves and discuss what happens, just as I did in problem 5.

dxdt6[x_,y_]:=(1-x/2-y/2)x;
dydt6[x_,y_]:=(2x-1)y/4;
slope6[x_,y_]:=((2x-1)y/4)/((1-x/2-y/2)x);

dt=0.01;
t0=0.0;
tEnd=20.0;

rk6[x0_,y0_]:=rungekutta2[dxdt6,dydt6,{t0,x0,y0},tEnd,dt];
curves6a={
rk6[3,3],rk6[2,1.5],rk6[0.1,0.1],
rk6[.1,.3],rk6[5,1],
rk6[3,.1],rk6[2,.1]
};

tEnd=1.0;
curves6b={
rk6[-.5,0],rk6[-.5,2],rk6[3,-.1],
rk6[-.1,-0.1],
rk6[0,-.5],rk6[-.5,-.5]
};

Show[
directionField[slope6, -1,5, -1,5,"x","y","problem 6"],
ListPlot[curves6a,PlotStyle->Red,Joined->True],
ListPlot[curves6b,PlotStyle->Green,Joined->True]
]

This phase space is more complicated than the others on the test.

I've plotted a number of curves above, in two colors; red ones starting with positive values, and green ones starting with x0 or y0 less than zero.  The green ones have a short ending time, to avoid blowing up.

Qualitatively, the curves that start in the positive quadrant all spiral in towards (0.5, 1.5), while the ones that start in other quadrants head out towards infinity.  

Analytically, it has 3 critical points : (0,0), (2,0), and (1/2,3/2),
all found by setting both derivatives to zero and solving the resulting equations;
in other words, curves that start at any of those three points stay motionless.

As shown below, two of those are unstable, one is stable.

Asking Mathematica to find these points can be done like this :

{crit1, crit2, crit3}=Solve[{dxdt6[x,y]==0,dydt6[x,y]==0},{x,y}]

Or that can be done by hand, by solving the two sets of equations simultaneously.

At each critical point , the equations can be linearized by approximating
the equations near each critical point.  The linear matrix is the change 
of dx/dt and dy/dt with respect to both x and y :

matrix6={
{D[dxdt6[x,y],x],D[dxdt6[x,y],y]},
{D[dydt6[x,y],x],D[dydt6[x,y],y]}}

{{1-x-(y/2),-(x/2)},{(y/2),(1/4) (-1+2 x)}}//MatrixForm

Here's the matrix at the (0,0) critical point, crit1 :

matrix6/.crit1//MatrixForm

Its eigens are

Eigenvalues[matrix6/.crit1]

So it's a saddle point.  This is consistent with the plot.

Likewise

Eigenvalues[matrix6/.crit2]

crit2

This one is an inward spiral, again consisten with the plot.  
Most of the phase points fall into this well, I think; 
it would take more exploration to know.

Eigenvalues[matrix6/.crit3]

The third critical point is another unstable saddle point.

And this is just what the text says in it's solutions at the end : 
 (0,0) is a saddle point, (2,0) is a saddle point, and (1/2, 3/2) is a stable inwards spiral.