Supplement for Lab 4: Ordinary Differential Equationsrestart;In this week's lab, you will learn to solve some basic systems of ordinary differential equations. These are the kinds of problems that computers can excel at because there are rote techniques for attacking the problems, but actual solutions often involve a lot of algebraic busywork.Some of the work you will be doing in this week's lab require loading some additional packages.with(plots):
<Text-field style="Heading 1" layout="Heading 1"></Text-field>How to Enter a Derivative in Maple
First of all, to enter differntial equations in Maple requires that you use the diff(y(t),t) differentiation function. This function takes the derivative of y(t) with respect to t. So we can write an expression for the differential equation ty' + 3y = 3ln(t) ast*diff(y(t),t)+3*y(t)=2*ln(t);We would often like to refer to this equation directly with a name, so I can instead use the Maple commandeqn1 := t*diff(y(t),t)+3*y(t)=2*ln(t);Notice that in the above equation, I have defined the entire differential equation ty' + 3y = 3ln(t) as eqn1 using the ":=" assignment operator versus the "=" equality operator. So what if I need higher order derivatives? Well, the diff function is meant to allow for any order of derivative, and even mixed partial derivatives, by simply adding the variables the derivative is taken against. So for example, for Newton's Second Law which has a second order derivative in time, I can enter diff(x(t),t,t) as the second order derivative of position with time, and so:NL2 := m*diff(x(t),t,t) = F(t);An alternate way to enter higher order derivatives is shown below:NL2_alternate := m*diff(x(t), t$2) - F(t);
For more help on diff, you can enter the following command (works for help on other commands too):?diff
<Text-field style="Heading 1" layout="Heading 1"></Text-field>Solving Differential Equations Analytically
As we have shown in class, some differential equations are relatively simple to attack and solve, such as the linear ordinary differential equations with constant coefficients. eqn2 := diff(y(t),t,t) + y(t) = cos(3*t);We can solve this differential equation using the dsolve function and telling it we want to solve for y(t) in eqn2:dsolve(eqn2,y(t));What does this result mean? It means that the solution of eqn2 has a sine and cosine components, with coefficients _C2 and _C1 respectively as well as an additional sine term with 3 times the frequency. _C1 and _C2 are essentially constants of integration and can only be solved by considering initial conditions on the system. Let's say we knew the system started from rest at the origin such that y(0) = 0 and y'(0) = 0. I can enter these initial conditions by changing the form of the dsolve function to one that allows us to also enter initial conditions. We will use the initial conditions y(0) = 0 and D(y)(0) = 0 where the second expression here is saying the first derivative of y at t=0 is 0.dsolve({eqn2,y(0) = 0, D(y)(0) = 0},y(t));The solution in this case is just the sum of two cosines with equal amplitude (opposite sign) but one is oscillating at triple the frequency of the other. We can plot this up by first defining the right hand side as a new function and then just plotting that function. In this case I will use the ditto operator % to grab the results from the last evaluated expression:eqn2_soln := rhs(%);plot(eqn2_soln,t=0..20);We can see more clearly what is going on here by plotting up the seperate cosine components of the solution on top of the entire solution:plot([eqn2_soln,1/8*cos(t),-1/8*cos(3*t)],t=0..20,color=[red,black,blue],legend=["Complete Solution","cos(t) Component","cos(3t) Component"]);I can see what effect varying the initial conditions has on the solution by starting off this solution with the object at the origin, but moving upward (in the process I will show another way to define initial conditions):ics := y(0) = 0, D(y)(0) = 0.5;dsolve({eqn2,ics},y(t));eqn2_2nd_soln := rhs(%);Notice the new cosine term that wasn't present in our previous "starting at rest" solution. I plot the two solutions against each other below so you can see the difference an initial velocity can make to the amplitude of the solution.plot([eqn2_soln,eqn2_2nd_soln],t=1..20,color=[red,blue],legend=["Start from Rest Soln","Swift Kick Solution"]);
Sometimes the results of running dsolve are not quite as pretty. We can try an equation that might look strangely familiar if you have been attending lecture:eqn3 := diff(x(t),t$5) + 17*diff(x(t),t$4) - 2*diff(x(t),t$3) +42*diff(x(t),t) - 10*x(t) = 0;
We can solve this differential equation using the dsolve function and we find:dsolve(eqn3);
Wow! That's ugly. What does it tell us? It tells us that the solution is the superposition of 5 exponential terms exp(at), each with a constant of integration _C_a where the five values of a are the five roots of z^5 +17*z^4 - 2*z^3 +42*z - 10 = 0. Which is about what we said in class! I can try expanding this out a bit by evaluating the series using the "simplification rules for algebraic numbers and functions"... in other words:
r2 := x(t) = rhs(%);r3 := evala(Simplify(value(r2)));r4 := evalf(r3);Now you know why Dr. Craig didn't do this on the board.... because the solution to this is "just some algebra" (Dr. Craig's words).
I can check if this is a solution I can use the odetest function which should evaluate to exactly zero if this is a solution. I can test the symbolic form with all the RootOf statements:odetest( r3, eqn3 );And I can test the evalf form which has been evaluated numerically and so some of the issues that we will be discussing shortly with computers rounding the results sets in.odetest( r4, eqn3 );Oops, looks like the rounding keeps the odetest function from returning exactly zero although some of those values are pretty small!
<Text-field style="Heading 1" layout="Heading 1"></Text-field>Solving Differential Equations Numerically
Here's another simple equation that will cause Maple to have a meltdown:eqn4 := diff(y(x),x,x) + diff(y(x),x) + x^5*y(x) = 0;dsolve(eqn4,y(x));You may have noticed Maple sat there and churned for a bit before returning essentially a restatement of the original problem. There is no exact analytical solution to this differential equation. It has to be tackled numerically. Again, Maple does a pretty good job with this. Let me set some initial conditions and this time I will numerically solve eqn4:ics4 := y(0) = 0, D(y)(0) = 3;eqn4_soln := dsolve({eqn4,ics4},y(x),numeric);The result looks a bit odd, its a Maple procedure, basically a string of commands, that can be called to give the numerical value of the solution at different points. For example, at x=4 the solution:eqn4_soln(4);have values of approximately y = -0.0905 and y' = 3.706.
Of course, in addition to seeing the value at a given point, we can plot the solution of y(x) as a function of x:odeplot(eqn4_soln,[x,y(x)],0..5);The plot looks a bit ratty, in part because the function varies rapidly with x and so the plot is just grabbing a few values of x, finding y(x), and "playing connect the dots." This can be a dangerous approach in that it can obsure the details. I can clean up the appearance of this plot immensely by increasing the number of points plotted:odeplot(eqn4_soln,[x,y(x)],0..5,numpoints=500);Wow, you can see a lot more of the details here.
I can also add the derivative to the same plot if I want using:odeplot(eqn4_soln,[[x,y(x)],[x,diff(y(x),x)]],0..5,numpoints=500,legend=["y","dy/dx"]);Does this plot make sense to you?
<Text-field style="Heading 1" layout="Heading 1"></Text-field>Animating a Plot by varying one parameter in the plot
In many differential equations, there is a parameter that determines the behavior of the solution, and Maple can be very helpful in examining those solutions. For example, consider the differential equation eqnvar := diff(x(t),t$2)+fb*x(t)=0;icsvar:= x(0)=3,D(x)(0)=0;The solution to this depends on the value of f. Let's find the solution:dsolve({eqnvar,icsvar},x(t));solnvar:=rhs(%);Now we will plot it as an animation, with different frames representing different values of fb. To see the animation, click on the plot and press the "play" button that appears in the main toolbar.animate(plot,[solnvar,t=0..10,y=-5..5],fb=-2..2);