<Text-field layout="Heading 1" style="Heading 1">Chapter 7: Ordinary differential equations</Text-field>restart: with(DEtools):
<Text-field layout="Heading 2" style="Heading 2">First order differential equations</Text-field>restart; with(DEtools):Radioactive decay:E1:= diff(N(t),t)= -g*N(t);General solution:dsolve( E1, N(t) );With initial condition NiMvLSUiTkc2IyIiISZGJUYm :dsolve({E1,N(0)=No},N(t));Before you can use the solution, you must capture it first.
<Text-field layout="Heading 3" style="Heading 3">Using rhs to get an expression</Text-field>S1:= dsolve( { E1, N(0)=No }, N(t) );Nsol:= rhs(S1);Plot:No:=5; g:=3; plot( Nsol, t=0..3);Note that N(t) is not defined:N(t);
<Text-field layout="Heading 3" style="Heading 3">Using rhs to get a function</Text-field>S1:= dsolve( { E1, N(0)=No }, N(t) );This doesn't work:Nsol:= t-> rhs(S1);Nsol(1); Nsol(2);This works as long as t is not assigned:Nsol:= x-> subs( t=x, rhs(S1) );Nsol(1); Nsol(2);Plot:No:=5; g:=3; plot( Nsol(t), t=0..3);Assigning t ruins Nsol:t:=3; Nsol(1); Nsol(2);Unassigning t gets its back:t:='t'; Nsol(1); Nsol(2);
<Text-field layout="Heading 3" style="Heading 3">Using assign</Text-field>N(t) is not yet defined:N(t);Now it is:assign( S1 ); N(t);Plot:plot(N(t),t=0..3);Differentiate, integrate:diff(N(t),t); int( N(t), t=0..3 );Note that N(t) is not exactly a function, so use it with care.These don't workN(5.); evalf(N(5.));This works:subs( t=5, N(t) ); evalf(%);Assigning t ruins N(t):t:=5; N(t); diff( N(t),t );Unassigning t cures it:t:='t'; N(t); diff( N(t),t );
<Text-field bookmark="Second order differential equation (top)" layout="Heading 2" style="Heading 2">Second order differential equations</Text-field>
<Text-field layout="Heading 3" style="Heading 3"> Newton's law of motion for a falling object</Text-field>restart: with(DEtools):E2:= diff( y(t),t\$2 ) =-g;Solution with initial conditions NiMvLSUieUc2IyIiISZGJUYm and NiMvLS0lJWRpZmZHNiQlInlHJSJ0RzYjIiIhJiUidkdGKg== :S2:= dsolve( { E2, y(0)=y, D(y)(0)=v }, y(t) );
<Text-field layout="Heading 3" style="Heading 3"> Harmonic oscillator</Text-field>restart: with(DEtools):E3:= diff(y(t),t\$2) = -omega^2*y(t);dsolve( E3, y(t) );Using copy & paste, we can re-write the solution as:ysol:= A*sin(omega*t)+B*cos(omega*t);To get a solution with initial conditions y(0)=1 and dy/dt(0)=2: s:= dsolve( { E3, y(0)=1, D(y)(0)=2 }, y(t) );assign(s); y(t);Plot:omega:=1; plot( y(t), t=0..20 );
<Text-field layout="Heading 3" style="Heading 3"> Damped harmonic oscillator</Text-field>restart:with(DEtools):Eq:= diff(y(t),t\$2) + diff(y(t),t)/tau + omega^2*y(t)=0;dsolve(Eq,y(t));
<Text-field layout="Heading 4" style="Heading 4"> Over-damped case</Text-field>NiMvJSJURyooIiIjIiIiJSNQaUdGJyUmb21lZ2FHISIi >> NiMlJHRhdUc= :omega:=2*Pi; tau:=.05;Y:= dsolve( { Eq, y(0)=1, D(y)(0)=0 },y(t) );plot( rhs(Y),t=0..20 );
<Text-field layout="Heading 4" style="Heading 4"> Under-damped case</Text-field>NiMvJSJURyooIiIjIiIiJSNQaUdGJyUmb21lZ2FHISIi << NiMlJHRhdUc= :omega:=2*Pi; tau:=20;Y:= dsolve( { Eq, y(0)=1, D(y)(0)=0 },y(t) );plot( rhs(Y),t=0..20 );
<Text-field layout="Heading 4" style="Heading 4"> Critical damping</Text-field>NiMvJSJURyooIiIjIiIiJSNQaUdGJyUmb21lZ2FHISIi = NiMqKCIiJSIiIiUjUGlHRiUlJHRhdUdGJQ== :omega:=2*Pi; tau:=1/(2*omega);Y:= dsolve( { Eq, y(0)=1, D(y)(0)=0 }, y(t) );plot( rhs(Y),t=0..20 );
<Text-field bookmark="Systems of differential equations (top)" layout="Heading 2" style="Heading 2">Systems of differential equations</Text-field>If we call the acceleration dv/dt and the velocity dx/dt, then the harmonic oscillator equation can be written as NiMvLSUlZGlmZkc2JC0lInZHNiMlInRHRiosJComJSZvbWVnYUciIiMtJSJ4R0YpIiIiISIi NiMvLSUlZGlmZkc2JC0lInhHNiMlInRHRiotJSJ2R0Yprestart: with(DEtools):Eq1:= diff(v(t),t) =-omega^2*x(t);Eq2:= diff(x(t),t) = v(t);XV:= dsolve( { Eq1, Eq2, x(0)=1, v(0)=0 }, {x(t),v(t)} );vsol:= rhs(XV);xsol:= rhs(XV);Plot:omega:=1; plot([xsol,vsol],t=0..20, color=[red,blue] );
<Text-field layout="Heading 2" style="Heading 2">Numerically solving differential equations</Text-field>F:= dsolve( ... type=numeric).
<Text-field layout="Heading 3" style="Heading 3"> Harmonic oscillator equations</Text-field>restart: omega:=1.;Eqx:= diff(x(t),t) = v(t);Eqv:= diff(v(t),t) =-omega^2*x(t);Initial conditionsIC:= x(0)=1.,v(0)=0.;Numerical solution:XV:= dsolve( { Eqx, Eqv, IC }, {x(t),v(t)}, type=numeric );Solution at t=2 :XV(2);Compare with the exact solution x(t)=cos(t):Digits:=18: cos(2.); Digits:=10:More accurate solution (use ?dsolve,numeric to see more):XV:= dsolve( {Eqx,Eqv,IC}, {x(t),v(t)}, type=numeric, method=rkf45, abserr=1e-10 );XV(2);Note: five or six significant figures is usually fine.
<Text-field layout="Heading 3" style="Heading 3"> odeplots </Text-field>To plot the solutions obtained in the last section, we need to use odeplot.First, the error tolerance of the solutions must be reduced:XV:= dsolve( {Eqx,Eqv,IC}, {x(t),v(t)}, type=numeric, method=rkf45, abserr=1e-6 );Then we plot:with(plots):odeplot( XV, [t,x(t)], 0..20, numpoints=500, color=navy );To plot both x(t) and v(t):odeplot( XV, [[t,x(t),color=red],[t,v(t),color=navy]], 0..20, numpoints=200 );
<Text-field layout="Heading 2" style="Heading 2">How does a differential equation mean?</Text-field>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. Explosive-growth:s:= dsolve( { diff(y(x),x)=y(x)^2, y(0)=1 },y(x) );plot( rhs(s), x=0..1 );
<Text-field layout="Heading 2" style="Heading 2">Nonlinear dynamics and chaos (advanced topic)</Text-field>restart; with(plots): with(DEtools):
<Text-field layout="Heading 3" style="Heading 3"> Lorentz equations </Text-field>Lorentz equations (for modelling the weather):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);If we choose NiMvJSZzaWdtYUciIzU= and NiMvJSJiRyomIiIpIiIiIiIkISIi , and then vary NiMlInJH from 0 up to about 30, a transition to chaos occurs.sigma:=10; b:=8/3;NiMvJSJyRyIjTg== :r:=35; s:= dsolve( {Eq1,Eq2,Eq3, x(0)=1,y(0)=1,z(0)=1}, {x(t),y(t),z(t)}, type=numeric );Plot:odeplot( s, [ [t,x(t),color=red], [t,y(t),color=blue], [t,z(t),color=tan] ], 0..20, numpoints=200 );NiMvJSJyRyIjNQ== :r:=10; s:= dsolve( {Eq1,Eq2,Eq3, x(0)=1,y(0)=1,z(0)=1}, {x(t),y(t),z(t)}, type=numeric );Plot:odeplot( s, [ [t,x(t),color=red], [t,y(t),color=blue], [t,z(t),color=tan] ], 0..20, numpoints=200 );NiMvJSJyRyIjRw== :r:=28; s:= dsolve( {Eq1,Eq2,Eq3, x(0)=1,y(0)=1,z(0)=1}, {x(t),y(t),z(t)}, type=numeric );Plot:odeplot( s, [ [t,x(t),color=red], [t,y(t),color=blue], [t,z(t),color=tan] ], 0..20, numpoints=200 );NiMvJSJyRyIjWA== :r:=45; s:= dsolve( {Eq1,Eq2,Eq3, x(0)=1,y(0)=1,z(0)=1}, {x(t),y(t),z(t)}, type=numeric );Plot:odeplot( s, [ [t,x(t),color=red], [t,y(t),color=blue], [t,z(t),color=tan] ], 0..20, numpoints=200 );NiMvJSJyRyIjWA== and NiMvLSUieEc2IyIiIS0lJkZsb2F0RzYkIiksKys1ISIo :r:=45; s:= dsolve( {Eq1,Eq2,Eq3, x(0)=1.0000001,y(0)=1,z(0)=1}, {x(t),y(t),z(t)}, type=numeric );Plot:odeplot( s, [ [t,x(t),color=red], [t,y(t),color=blue], [t,z(t),color=tan] ], 0..20, numpoints=200 );
<Text-field layout="Heading 3" style="Heading 3"> Spacecurves</Text-field>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 rr:=50;for i from 1 to numpoints do s:= dsolve( { Eq1,Eq2,Eq3, x(0)=r0, y(0)=r0, z(0)=r0}, {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 againod:Turn rf into a sequence for use with spacecurveL:= [ seq( rf[i], i=1..numpoints ) ]:Plot the trajectoryspacecurve( 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 characterfor 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 );