Chapter7.mws

• First order differential equations
• Second order differential equations
• Systems of differential equations
• Numerically solving differential equations
• How does a differential equation mean?
• Nonlinear dynamics and chaos (advanced topic)

---

Chapter 7: Ordinary differential equations

Differential equations are the language of physics and Maple can help you solve them. It is a good idea before starting a problem with a differential equation in it to load Maple's full set of differential equation tools, like this

> restart:with(DEtools):

>

First order differential equations

Debugging

> restart;with(DEtools):

Let's begin by asking Maple to solve the first differential equation you probably ever saw, the equation for radioactive decay. If the decay rate of the atoms in a sample is (g has units of 1/second), then the equation determining the number of atoms still left in the sample is

(Differential equations are built like ordinary equations in Maple, but using the diff command:)

> E1:=diff(N(t),t)=-g*N(t);

along with the initial condition that at time there were atoms in the sample. The problem is to find the function . Maple solves it this way using the command dsolve .

> dsolve(E1,N(t));

As expected, the solution is a dying exponential. Maple gives you an unknown constant out front because we didn't tell it what the initial conditions are yet. The constant is not too hard to find because at the exponential function is 1, so the constant must be . But dsolve will accept initial conditions as well and do even this little bit of work for you. You simply replace the first argument of dsolve with a set containing E1 and the initial condition, like this. (See ?dsolve,initial_conditions for more details)

> dsolve({E1,N(0)=No},N(t));

Ok, now suppose you wanted to get a plot of this solution. First you need to tell Maple what the values of and are like this

> No := 5;g:=3;

and then you do something that seems a little odd. You have to somehow grab the expression on the right-hand side of the solution from dsolve so you can pass it into plot. One way to do this would be to start the plot command like this,

plot(

then stop, use the mouse to highlight the expression on the right-hand side of N(t)=... above and paste it in. But this seems awkward; surely Maple will give you access to the solution it just found? It will, and here is one way that works. Instead of just solving the differential equation to the screen, assign the solution to a variable, like this

> S1:=dsolve({E1,N(0)=No},N(t));

You can now refer to this solution by its name, S1, and in particular you can refer to the solution itself by rhs(S1) , like this

> plot(rhs(S1),t=0..3);

Note that the command rhs picks off the right-hand side of an equation. You can probably guess what lhs does.

>

Another way to plot this result is to tell Maple to invent a variable with the name N(t) which can be used as an expression in the plot command. The syntax for this is simply assign(S1); .

> assign(S1);

and then we can use

> plot(N(t),t=0..3);

The question naturally arises, "What exactly has become of N?" If you try to look at it this way

> N;

there isn't anything there. The problem is that Maple did something perfectly natural for it, but rather odd for us. It took the result for S1, which is quite literally and "sort of " made an expression called N(t). Look at

> N(t);

It makes sense to differentiate this thing

> diff(N(t),t);

It makes sense to solve for when it is 1

> solve(N(t)=1,t);

But you can easily see that it isn't a function that you can plug numbers into by doing this:

> N(5);

And if you try to give t a value and then evaluate the expression N(t), not only does it not work, but it ruins N(t) as well:

> t:=3;

> evalf(N(t));

> solve(N(t)=1,t);

Error, (in solve) a constant is invalid as a variable, 3

The problem seems to be that the symbol N(t) that comes out of the assign is a funny mix of expression and function, so if you want to use this assign syntax instead of rhs , and it often does seem to be more natural to do so, be very careful about assigning values to the things that appear in the expression. And also notice that you can't go back to dsolve again with different initial conditions either, because N(t) has been ruined--it's time to restart .

> S1:=dsolve({E1,N(0)=3},N(t));

Error, wrong number (or type) of parameters in function diff

>

This may seem so troublesome that you might never want to us assign , but assign makes the plot and other commands look so much cleaner than they do with rhs that I think it is still worth it sometimes. But don't plan on using it exclusively; rhs is often the safer thing to use.

One final note about this: in the plot syntax, in the solve syntax, and in the differentiation syntax when it looks like we are using a Maple function N and invoking it with N(t) we are actually using a Maple expression (sort of) called N(t). If this is confusing (and it certainly confuses me) it might help to review the difference between expressions and functions way back in Chapter 1. Expressions and functions

Problem 7.1

The differential equation governing the current in a circuit containing a battery of emf , resistance , and inductance is

.

Solve this differential equation for the current and plot the result if the initial current is zero, H, , and V.

----------------------------------------------------------------------------

Problem 7.2

Here is a problem involving the idea of "explosive growth". Suppose that you have discovered some process in which the rate of growth of the quantity is not proportional to itself, as in exponential growth, but is instead proportional to some power of , say or . Use Maple to explore these two cases by having it solve the differential equation

for and and with the initial condition . Then plot the solution over a big enough time interval that you can see the explosion.

----------------------------------------------------------------------------

Problem 7.3

Consider the relatively innocuous first-order differential equation

with .

Ask Maple to solve it and see what happens. In a later section on solving differential equations numerically we will try to solve this equation again.

Go to top of section

Second order differential equations

Debugging

Second order equations are very common in physics because Newton's second law is second order and because diffusion, waves, and quantum mechanics all have governing equations containing second derivatives. Let's start with Newton's law of motion, and just as a warm up, let's do a falling object. The differential equation is just

> restart:with(DEtools):

> E2:=diff(y(t),t$2) = -g;

> S2:=dsolve({E2,y(0)=y[0],D(y)(0)=v[0]},y(t));

Well, that should be a familiar result. Notice a few things about the way this was done. To enter an initial condition on the derivative of y use D(y) , which means dy/dt ( D is Maple's generic derivative symbol) then tack on (0) to tell it that the derivative is evaluated at t=0. Yes, this means that diff and D both mean derivative, and there is an important and subtle difference between them that is hard to remember; it is discussed in the section on derivatives in Chapter 3( diff and D). So either memorize the syntax of the dsolve command above, or remember where it is in this document so you can easily find it when you need to solve problems like this.

Here's another famous second order differential equation, the harmonic oscillator.

> E3:=diff(y(t),t$2)=-omega^2*y(t);

Maple solves it in its generic form with arbitrary constants this way

> dsolve(E3,y(t));

Warning: this form of the solution is dangerous because there is no guarantee that every time you run a worksheet with this dsolve command in it that _C1 will go with sine and _C2 will go with cosine. Maple has its own internal logic about how to make this choice, but it can turn out to look pretty random in more complicated differential equations than this one. So don't let Maple make the choice. Once you get this form of the solution, use the mouse to copy the solution onto a new line and choose your own unknown coefficients like this (and get rid of the awkward y(t) = and replace it with a genuine assignment as well:)

> ysol:= A*sin(omega*t)+B*cos(omega*t);

Any calculations you do to determine A and B will now work the same way every time because you have determined for all time that A goes with sine and B goes with cosine.

If the constants are to be determined by initial conditions there is a form of the dsolve command that will determine them automatically. For instance, suppose that you want to solve the differential equation E3 with initial conditions y(0)=1 and dy/dt(0)=2. Then you can use

> s:=dsolve({E3,y(0)=1,D(y)(0)=2},y(t));

Now we have Maple give us the expression y(t) by using assign

> assign(s);

And then we can plot the solution

> plot(y(t),t=0..20);

Hmmm; that didn't work too well. The plot frame is there, but no function. And, unfortunately, no hint about what went wrong. To debug this problem it would be nice to put some values for t into y(t), but remember that y(t) looks like an expression, but isn't quite, and that assigning t a value will ruin it. So let's just look at y(t) this way and see if we can see what's wrong

> y(t);

The plot command will supply numerical values for t, but the fact that is just sitting in this expression without having a numerical value tells us what is wrong; we haven't given a numerical value! Ok, so now do this

> omega:=1;plot(y(t),t=0..20);

Now this problem was simple enough that you could have easily done it yourself once you had taken a class about differential equations. But Maple can do harder problems with ease as well. Let's add some damping to the harmonic oscillator by putting in a damping force of the form , where is the characteristic time for damping to occur. The new equation of motion is now

> restart:with(DEtools):

> Eq:=diff(y(t),t$2) + diff(y(t),t)/tau + omega^2*y(t)=0;

and we get the general solution this way

> dsolve(Eq,y(t));

>

Well, that looks impressive. But if it's supposed to be an oscillator that damps, where are the sines and cosines? The problem is that Maple doesn't yet know how big and are. Let's think physically for a minute. Suppose we have a natural frequency of about , so the period is 1 second. Now let's suppose that we have this oscillator (think of a pendulum) in motor oil at 50 degrees below zero so that the characteristic damping time is 0.05 seconds (so much friction that motion stops in 0.05 seconds). Your intuition should tell you that if you pull the pendulum back and release it, you are not going to see any swinging; all that is going to happen is that the pendulum will slowly ooze its way into the vertical position and stay there. Let's see what Maple says

> omega:=2*Pi;tau:=.05;

> Y:=dsolve({Eq,y(0)=1,D(y)(0)=0},y(t));

> plot(rhs(Y),t=0..20);

Ok, Maple agrees with your intuition.

>

Now let's change the damping and imagine what will happen. Your intuition should tell you that if you warm the oil up, or use WD-40 instead, or maybe even just air, the pendulum will start to swing and only damp slowly. For example, if then Maple gives us

> omega:=2*Pi;tau:=20;

> Y:=dsolve({Eq,y(0)=1,D(y)(0)=0},y(t));

> plot(rhs(Y),t=0..20);

Well, now it swings and only slowly damps. By what miracle did the exponential functions in the original solution above become sines and cosines? Remember Euler's formula: ; in this oscillating case the constants in the exponential become imaginary, and that's how the formula changes from decay to wiggles.

>

This brings up the question of what value of makes the transition from pure damping to damped oscillation. Look back at the general form and you will see that there is a square root

in the exponentials. If the argument of this square root is positive, then both of the fundamental solutions in the sum are decaying exponentials and we have damping. If the argument of the square root is negative, then the square root is imaginary, the exponentials are imaginary exponentials, i.e., sines and cosines, and we have oscillation. The transition between the two is when the argument of the square roots is zero, i.e., when . This special case is called critical damping.

Problem 7.4

Have Maple solve the damped harmonic oscillator differential equation for the case of critical damping with and plot the function y(t) from t=0 to t=2 with initial conditions y(0)=1 and dy/dt(0)=0. The plot will look just like damping with no hint that you are right on the border between damping and oscillation. To verify that you are right on the border increase by 2% and make the plot again. Window the graph down enough that you can see that y(t) goes a little bit negative, indicating that oscillation has begun. You will have to do some serious windowing to see this effect.

----------------------------------------------------------------------------

Problem 7. 5

One last harmonic oscillator problem:

(a) Let's add a sinusoidal driving force at a frequency a little bit different from the natural frequency , say . Have Maple solve this differential equation and plot it from t=0 to t=200 with and . Use completely-at-rest initial conditions, y(0)=0 and dy/dt(0)=0.

The plot should show the classic pattern of the driven-damped harmonic oscillator: an initial period of start-up with some beating between the natural and driving frequencies, then a transition to a state of oscillation with y(t) matched to the driving frequency. Verify graphically that the final oscillation frequency of y(t) is .

(b) Study resonance in this system by making plots of vs. in the following cases:

(i) : , , and .

(ii) : , , and .

(iii) : , , and .

Verify that as the damping time gets longer, the resonant amplitude gets larger and the resonance gets narrower.

(c) Find formulas for the steady state motion of this driven-damped oscillator by substituting the guess

into the differential equation in part (a). Then require the and in the resulting equation to be separately equal (they are independent functions) to get two equations for the two unknowns and . Solve these two equations and verify that the formula for agrees with the steady state amplitudes you found in part (b). Finally, make plots of and for the three choices of used in part (b) and verify again that a long damping time leads to larger and sharper resonance.

----------------------------------------------------------------------------

Problem 7.6

Use Maple to solve the following second-order differential equations in general form.

(a)

(b)

(c)

Go to top of section

Systems of differential equations

Debugging

It is often more convenient to deal with systems of differential equations than with second, third, or higher order differential equations. Indeed, many numerical methods require that you write your differential equation as a system of first order differential equations. Here's a simple example. If we call the acceleration dv/dt and the velocity dx/dt, then the harmonic oscillator equation can be written as the following first order set

The Maple command dsolve handles systems like this by asking for a set {sets have curly brackets} of equations and initial conditions, a set of functions to be solved for, and perhaps extra arguments, like this

> restart:with(DEtools):

> Eq1:=diff(v(t),t)=-omega^2*x(t);Eq2:=diff(x(t),t)=v(t);

> dsolve({Eq1,Eq2,x(0)=1,v(0)=0},{x(t),v(t)});

Notice that this is actually a nice form to have the answer in because we get both the position and the velocity simultaneously. And another good thing is that you don't have to remember the D notation for derivative initial conditions. To plot x(t) and v(t) we first need to give a value, then solve again

> omega:=1;XV:=dsolve({Eq1,Eq2,x(0)=1,v(0)=0},{x(t),v(t)});

then make the plots like this using the assign command again

> assign(XV);

> plot([x(t),v(t)],t=0..20);

>

Problem 7.7

Here are the equations of motion for a charged particle of mass m and charge q moving in the field of a constant magnetic field B in the z-direction and a constant electric field E in the y-direction which come from Newton's second law using the Lorentz force law: . (Note: the phrase equation of motion has special meaning for physicists and it takes a while for beginning physicists to figure out what we are talking about. It seems natural at first to think that the phrase "equation of motion" must mean the equations that describe the motion, like , or . This is not what physicists mean when they say "equation of motion". By equation of motion we mean the differential equation which, when solved, gives the formulas that describe the particle motion. If you can get this straight early on it will save you a lot of trouble.)

Define these equations in Maple and have Maple solve them with q=-1.6e-19, m=9.11e-31, B=0.2, E=3000 with initial conditions

, , , . Make plots in the x-y plane of the particle orbits and (they must be parametric plots with t as the parameter). You will need to choose a time interval for plotting and if you just choose something simple like 1 or 2, you will get garbage. Think about the physics of the problem; in particular, think about the cyclotron frequency: . Also make x-y plane plots of and with the scaling=constrained option. In this way you will be able to see that the velocity goes around in a circle (cyclotron motion) but with an offset. Verify that the velocity offset (called the "drift velocity") is given by .

Go to top of section

Numerically solving differential equations

Debugging

Most of the time Maple will not be able to solve your differential equation for you in terms of functions that it knows. When this happens you will want to solve the equation numerically, i.e., perhaps make a plot of the solution, find out what the solution is at x=2, etc.. You tell Maple to solve the differential equation numerically by using dsolve with the option type=numeric , as shown in the example below. But before we start doing this let me warn you. (1) "Numerical" implies that floating point numbers (numbers with decimal points) are used. If you give Maple integers for values in initial conditions or constant definitions, dsolve(....type=numeric) may fail. (2) Since the whole point of doing a numerical calculation is to get a number and since Maple is loath to ever give you numbers, doing useful things with what Maple returns as the numerical solution of a differential equation is a bit awkward. You are likely to become more frustrated in this section of the text than any other so far. (3) You will have much better luck at getting Maple to give you what you want if you define your differential equations as first order sets than if you use 2nd, and higher, order derivatives. (4) Assign the result of dsolve to a variable name so you can refer to it later, i.e., use F:=dsolve(...type=numeric).

OK, let's start by doing the harmonic oscillator numerically. First define the equations for x(t) and v(t).

> restart:omega:=1.;

> Eqx:=diff(x(t),t)=v(t);Eqv:=diff(v(t),t)=-omega^2*x(t);

then define the initial conditions

> IC:=x(0)=1.,v(0)=0.;

and finally invoke dsolve with the numerical option set to get the numerical solution

> XV:=dsolve({Eqx,Eqv,IC},{x(t),v(t)},type=numeric);

XV is now a procedure that expects to be given a value of t so that it can return an approximate solution to x(t) and v(t), like this

> XV(2.);

Notice that the answer came out apparently accurate to 18 decimal places. In this case we know the answer is just , so we can check it (after increasing the number of digits Maple keeps in floating point calculations).

> Digits:=18:cos(2.);

Apparently Maple lied--the result of dsolve is pretty accurate (7 decimal places) but not as accurate as the numbers it returned would indicate.

>

We can make it more accurate at the expense of having the solution take longer by specifying the numerical method to use (use ?dsolve,numeric to see what the choices are; rkf45 is a pretty good choice; if you have trouble, try another method) and setting the allowable error like this

> XV:=dsolve({Eqx,Eqv,IC},{x(t),v(t)},type=numeric,method=rkf45,abserr=1e-10);

> XV(2);

>

So this is better, but Maple won't let you improve it forever. Try changing the error to 1e-15 in the command above and see what happens. (Don't leave it this small in what follows or you will work your computer to death.)

And now a word of advice: don't try to compute to an accuracy of 1 part in 1000000000000. Five or six significant figures is usually fine, although it is a good idea to do a numerical calculation twice with different error tolerances to make sure the setting you are using most of the time is OK for the problem at hand.

All right, Maple has given us a procedure called XV that solves this differential equation numerically. Now suppose we want to plot it. Here's one thing that won't work: you can't use assign(XV) like we did with symbolic solutions. The easiest way that does work is to use the Maple command odeplots which expects to be given output from dsolve(....type=numeric) . Note that you have to use the command with(plots) to make this command available.

Here's another problem before we proceed. The error tolerance in the dsolve command above is set to 1e-10, which is so small that asking odeplot to find lots of points at this level is asking to be bored in front of a busy screen for a long time. So before we proceed, let's set the error to a more reasonable value for plotting.

> XV:=dsolve({Eqx,Eqv,IC},{x(t),v(t)},type=numeric,method=rkf45,abserr=1e-6);

In this case it would be used like this to plot x(t)

> with(plots):odeplot(XV,[t,x(t)],0..20,numpoints=500,color=navy);

Warning, the name changecoords has been redefined

or this way to plot both x(t) and v(t)

> odeplot(XV,[[t,x(t)],[t,v(t)]],0..20,numpoints=200);

Note two things. (1) odeplot takes all the usual plot options so you can make it look good and (2) if you have your error set really low and numpoints set really high, you are going to be waiting a long time for your plot to come out.

>

Problem 7.8

Here is problem 7.3 again, which Maple couldn't do symbolically. Do it numerically and plot y(x) from x=0 to x=20.

with .

----------------------------------------------------------------------------

Problem 7.9

Here is a real physics problem. In Physics 121 you did the problem of a hard-hit baseball, but because you did it without air friction you were playing baseball on the moon. Let's play ball in a real atmosphere now. The air-friction drag on a baseball is approximately given by the following formula

where is the drag coefficient (about 0.35 for a baseball), where is the density of air (about 1.2 kg/ ), where is the radius of the ball (about 0.037 m) and where is the vector velocity of the ball. The absolute value in this formula pretty much guarantees that Maple won't be giving you a formula for the solution of this problem. If we take the vertical direction to be y and the horizontal direction to be x then Newton's second law in vector-component form gives us the following first-order set of equations

where is the mass of the baseball (0.145 kg) and where is the acceleration of gravity (9.8 m/ ). Give this problem to Maple by coding the four differential equations and hit a home run with it using initial conditions

where m/s. Then use dsolve(...type=numeric) to solve it and make overlaid plots: (i) without air drag (set ) and (ii) with air drag, as discussed above. Plot y(t) vs. x(t) so you can see the path of the ball by using [x(t),y(t)] in odeplot instead of [t,y(t)]. Can you see that the trajectory with air friction is not a parabola? Power hitters in baseball say they would much rather play in Coors Field in Denver than in sea-level stadiums because it is so much easier to hit home runs. Do they know what they are talking about? The air pressure in Denver is about 10% lower than it is at sea level.

----------------------------------------------------------------------------

Advanced topic: shooting

Problem 7.10

Sometimes we don't want to just solve a differential equation with known initial values. Perhaps there is a parameter in the differential equation that is unknown and we need to determine its value based on the solution of the differential equation at some other point. Or perhaps we want to start out with part of the initial conditions given and then adjust the undetermined initial conditions until that the solution satisfies some other conditions at the far end. Adusting things so that you start at point a and then hit a target at point b sounds so much like artillery that techniques that solve this kind of problem are usually called shooting methods. When a symbolic solution can be found this is not so hard, but when the differential equation must be solved numerically it can be quite difficult. In fact, I can't make Maple do it with its built-in commands, so I will have to give you a procedure based on one written by Dr. Robert Israel to show you how such problems are solved. But first, let's choose a specific problem to solve.

Schrodinger's equation for the bound states of a particle in a harmonic oscillator potential well is this:

To make it a nice dimensionless problem we change the length variable from to by doing this:

where

and by changing into the eigenvalue by doing this:

.

With these changes of variable Schrodinger's equation becomes

.

The solutions of this equation come in two types: even functions in and odd functions in . Hence we will have two different kinds of initial conditions at :

even: and odd: and

(You may be wondering why we use 1 for the initial value or slope; why not .7, 1e-16, or something else? What we use doesn't matter. If you look at the differential equation every term has in it, so if you double some solution that you have found, the doubled will also be a solution, just twice as big. So since the size doesn't matter we just use 1. If you were going to do real physics, however, then at the end you would have to choose an overall multiplier for so that the total probability is 1. You will learn all about this in quantum mechanics.)

The only solutions for the wave function that make physical sense are the ones that go to zero as approaches . You will shortly see that this is only true for very special values of the parameter , so we need a procedure that does the following: (i) It must be able to solve the differential equation numerically and (ii) it must be able to adjust the parameter to make the wavefunction be small at large .

Well, here is the Maple code containing the shooting procedure, all set up to solve a simpler problem, which is similar to our quantum problem. The practice differential equation is (harmonic oscillator)

with and .

We want to adjust so that we meet the condition . Follow along and run the procedure below, then modify it to find the values of in Schrodinger's equation for the first 6 bound states (3 even and 3 odd). You will need to decide how far out to go in to approximate infinity. Don't overdo it! Somewhere in the 5-10 range is probably enough.

Start fresh and pull in odeplot

> restart:with(plots):

Warning, the name changecoords has been redefined

Define the differential equation and assign it to deqn. NOTE: the parameter you are going to vary to satisfy the shooting condition on the right end must be called lambda for this section of code to work.

> deqn:=diff(y(x),x$2)+lambda^2*y(x)=0;

Define the left and right endpoints of the interval

> a:=0;b:=1;

Define the initial conditions at the beginning of the interval.

> bc:=y(a)=0,D(y)(a)=1;

Set the value that you are shooting for, i.e, set the variable yright to the value that y(b) is supposed to be when you have the correct value of the parameter

> yright:=0;

Make a guess at the parameter to be varied and plot the numerical solution of the ode from a to b; Adjust the value of until you are close to the solution you want. (In the sample problem the parameter values that work are , , , etc..

> lambda:=2.5;

Note that the startinit=true option tells dsolve to always start with the specified initial conditions. This is important to have set when we use fsolve to find the value of param that gives the solution with y(b)=yright.

> S:=dsolve({deqn,bc},y(x),type=numeric,startinit=true);

> odeplot(S,[x,y(x)],a..b);

Go back up and try another value of until you get close to zero on the right.

I am now going to show you a short section of code that makes it possible to pass the result of dsolve into fsolve to let Maple find the correct value of . The way to make this work was shown to me by Dr. Robert Israel of the University of British Columbia.

>

Note that this code won't work unless the dummy argument of the procedure F given below (z in this case) is different from all of the variables in the differential equation, i.e., in our problem you can't use y or x in place of z below or you will get in trouble. And because has already been defined, you can't use it either; that's why I use the dummy variable z in place of . (But why would you want to use x or y, you ask? In this problem you wouldn't, but in a different problem you might be tempted, and I want you to know that it's a mortal sin.)

> F:=proc(z)

> local a,b,deqn,bc,S;

> deqn:=diff(y(x),x$2)+z^2*y(x)=0;

> a:=0;b:=1;

> bc:=y(a)=0,D(y)(a)=1;

> if type(z,numeric) then

> S:=dsolve({deqn,bc},y(x),type=numeric,startinit=true,abserr=1e-6);

> subs(S(b),y(x))

> else

> 'F(z)';

> end if;

> end;

F( ) now can be passed into fsolve to find a value of that satisfies F(b)=yright. In the paragraph below you will see how to do this.

Set Digits to a lower value so fsolve doesn't kill itself trying to get an answer correct to many significant figures when dsolve isn't accurate enough to deliver this much accuracy. If you don't get the right balance between abserr in the dsolve in the procedure and the Digits setting here, fsolve will run forever.

> Digits:=5;

> F(3.);

> lambda:=fsolve(F(z)=yright,z=3.0..3.2);

Using this value of , make a plot of the solution.

> S:=dsolve({deqn,bc},y(x),type=numeric):

> odeplot(S,[x,y(x)],a..b);

>

----------------------------------------------------------------------------------

Problem 7.11

The equation of motion of a relativistic harmonic oscillator is , where is the relativistic momentum , where is the particle velocity . (a) Write this equation of motion in the form of two coupled first order differential equations. The first one is easy because it has been given in the first line of this problem. The second one is of the form . Solve for in terms of to get this second equation. (b) Use dsolve(...type=numeric) to make a procedure that you can use in odeplot to make graphs of for the cases , .4, 1., 5. and . Verify graphically in the last case that the reason that x(t) looks so strange is that it is obeying the relativistic speed limit. In this last case make a plot of the particle speed u(t).

Go to top of section

How does a differential equation mean?

Debugging

It is important in physics to have an intuitive feel for what differential equations mean rather than just memorizing a bunch of techniques and formulas. Without this feel you won't be able to propose and refine mathematical models for physical processes. So let's look at a simple one and try to translate it into words:

Since is the slope of the function y(x) this equation says that the bigger y gets the bigger its slope gets. Let's try the two possible cases. (i) y(0) is positive. The differential equation then says that the slope is positive, so y is increasing. But if y increases its slope increases, making y increase more, making its slope increase more, etc.. So it is no surprise that the solution of this equation is a function like that gets huge as x increases. (ii) y(0) is negative. Now the differential equation says that the slope is negative, so y will have to decrease, i.e., become more negative than it was at x=0. But if y is more negative then the slope is more negative, making y even more negative, etc.. Again, it's no surprise that the solution is a strongly decreasing function like .

With this kind of qualitative analysis in mind, let's go back to an explosive-growth differential equation like the ones you did in the section on first order equations, say

.

Keeping in mind that with a first power on the right we got the exponential function, this equation says that if y starts out positive, y should increase even more than it did before, i.e., get bigger faster than the exponential function. That would have to be pretty impressive, and it is: it goes to infinity before x gets to infinity:

> s:=dsolve({diff(y(x),x)=y(x)^2,y(0)=1},y(x));

> plot(rhs(s),x=0..1);

Problem 7.12

Now think about what should happen if y(0) is negative. Make some rough sketches, then check your answer by changing the initial conditions in the example above.

----------------------------------------------------------------------------

>

You can play this game with second-order differential equations too. Let translate the simple harmonic oscillator equation

into words. We need to remember that the second derivative means the curvature of the function: a positive second derivative means that the function curves like the smiley face of someone who is always positive, while negative curvature means that it curves like a frowny face (easy to remember, no?) And if the second derivative is large in magnitude then the smile or frown is very narrow (like a piece of string suspended between its two ends from fingers held close together) and if the second derivative is small in magnitude it is like taking the same piece of string and stretching your arms apart to make a wide smile or frown. Ok, what does this equation say if y=1 to start? Well, the second derivative is negative, so the function y curves downward, making y smaller, which makes the frowniness smaller, but still negative, so y keeps curving downward until it crosses y=0. Then with y negative the differential equation says that the curvature is positive, making y start to smile and curve upward. It doesn't curve much at first because y is pretty small in magnitude, but eventually y will have a large enough negative value that y(x) turns into a full-fledged smile, stops going negative, and heads back up toward y=0 again. When it gets there y becomes positive, the function gets frowny and turns back around toward y=0, etc.. So it is no surprise that the solution of this equation is an oscillation, .

Problem 7.13

(a) Use this qualitative analysis to guess what the solution of this equation would look like, then check your answer with Maple.

with and .

(b) Use this qualitative analysis to guess what the solution of this equation would look like, then check your answer with Maple (you will need to use dsolve(....type=numeric) and odeplot ).

with and

(c) Same thing for this one

with and

----------------------------------------------------------------------------

Problem 7.14

Think again about the baseball problem, Problem 7.9. If you look at your trajectory, and expecially if you increase the initial speed of the bat up to something like 80 m/s, you will see that instead of a parabola you get a curve that starts like a parabola, but then drops almost vertically at the end. Use qualititative analysis on the x and y components of Newton's second law to explain why this occurs.

----------------------------------------------------------------------------

Problem 7.15

Consider again the harmonic-oscillator bound-state problem, Problem 7.10:

Think about what its solutions qualitatively look like both for small and for large and show that this qualitative analysis gives roughly the right answer for what the bound state eigenfunctions look like.

Go to top of section

Nonlinear dynamics and chaos (advanced topic)

Debugging

The use of computers over the past 40 years has completely revived the stodgy old subject of classical mechanics, and the system of equations that started the whole thing are these three coupled ode's discovered by Ed Lorentz while trying to model the weather.

> restart;with(plots):

> Eq1:=diff(x(t),t) = sigma*(y(t)-x(t));

> Eq2:=diff(y(t),t) = r*x(t)-y(t)-x(t)*z(t);

> Eq3:=diff(z(t),t) = x(t)*y(t)-b*z(t);

Warning, the name changecoords has been redefined

where the constants , , and can chosen however you like. For most choices of these parameters the solutions are not particularly interesting, but Lorentz discovered that there is a range of choices for which the solution behaves in a very strange way. If we choose

> sigma:=10;b:=8/3;

and then vary from 0 up to about 30, a transition occurs from ordinary boring behavior to something so strange that it has revolutionized many fields of science, including physics. To numerically solve this set for a specific choice of r we use

> r:=35;s:=dsolve({Eq1,Eq2,Eq3,x(0)=1,y(0)=1,z(0)=1},{x(t),y(t),z(t)},type=numeric);

and to see what the solutions look like we use odeplot

> odeplot(s,[[t,x(t)],[t,y(t)],[t,z(t)]],0..20,numpoints=200);

Take a few minutes now and run this execution group for the following values of : 10, 20, 28, 35, 45. For r above about 25 the oscillations form a never-repeating pattern (which you can sort of see if you run for a longer time in the odeplot ). And what's worse, if you use slightly different initial conditions you get a completely different pattern (this is the famous butterfly effect. To observe it change x(0)=1 in the initial conditions to x(0)=1.0000001.) This is completely unexpected. The equations that govern this behavior are completely deterministic and quite simple--how can random garbage come out of something this simple? It's truly amazing, and this kind of behavior is called dynamical chaos .

>

Here is another way to look at it. Instead of plotting x(t), let's plot the trajectory of the solution curve [x(t),y(t),z(t)] as a 3-dimensional curve using Maple's spacecurve plot command. To use it we need to make a loop that repeatedly calls dsolve to generate a sequence of 3-dimensional points along the system trajectory. This sequence of points is then passed into spacecurve , like this:

Put the initial values of (x,y,z) in r0, the starting time in t1, the ending time in t2 and the number of points along the trajectory in numpoints ; then calculate the time step between successive points.

> r0:=[1,1,1];t1:=0;t2:=10;numpoints:=400;dt:=(t2-t1)/numpoints;

Set the system parameter r

> r:=50;

>

Step the system along with repeated calls to dsolve , each one having for its initial conditions the ending conditions from the previous step. This is done with a for-do-end do loop, which we will study in Chapter 8. Note: this takes a while, so be patient.

> for i from 1 to numpoints do

> s:=dsolve({Eq1,Eq2,Eq3,x(0)=r0[1],y(0)=r0[2],z(0)=r0[3]},{x(t),y(t),z(t)},type=numeric):

> t0:=t1+(i-1)*dt:

> tf:=t0+dt:

> sol:=s(tf-t0):

Store (x,y,z) in rf[i] after each step

> rf[i]:=[rhs(op(2,sol)),rhs(op(3,sol)),rhs(op(4,sol))]:

Get the next initial condition

> r0:=rf[i]:

Bottom of the loop; go back up and do it again

> od:

Turn rf into a sequence for use with spacecurve

> L:=[seq(rf[i],i=1..numpoints)]:

Plot the trajec tory

> spacecurve(L,orientation=[45,60]);

>

Now make 100 spacecurve plots, each one with a slightly different viewing angle. Then we will animate the sequence so that the trajectory rotates on the screen, making it easy to see its 3-dimensional character

> for i from 1 to 100 do

> theta :=45:phi:=i*3.6:

>

> p[i]:=spacecurve(L,orientation=[theta,phi]):

>

> od:

Plot the animation. Note: when the picture comes up and the clock goes away, click on the picture to make the animation controls appear on the toolbar.

> plots[display]([seq(p[i],i=1..100)],insequence=true);

Of course you can just click on a single space curve plot and use the cursor to rotate the picture any way you want, but I wanted you to see how to animate a spacecurve.

>

Problem 7.16

Look at how the system trajectory changes character as it changes from going to a fixed point to wandering around on the so-called strange attractor displayed by spacecurve. The transition occurs as r increases past a number around 25, or so, although if you let the system run for a long time the attractor starts to show up at lower values of r. If you try to keep numpoints too high this will take forever, and if you don't run long enough you won't be able to see where the transition occurs; use numpoints=200 and t2=20.

Go to top of section

Go to top of chapter

>

>