Mathematica: Differential Equations
PHYS 460/660 @ UDelaware

Reference: Heikki Ruskeepaa, Mathematica Navigator (Academic Press, San Diego, 2004).

To open a section, double click on a cell bracket at the right side of the screen that has a hook on its bottom end like this ⥡

0. Initialization

ClearAll["Global`*"]

1. Symbolic Solution

({{Symbolic Solution for  Single ODE}, {({{Give the general solution, sol=y[t]/.Dsolve[eqn,y[t],t]}, {Solve an initial value problem, sol=y[t]/.Dsolve[eqn,y[a]==α,y[t],t]}, {Plot the solution of an initial value problem, Plot[sol,{t,a,b}]}})}})

Solve first-order ODE of the linear harmonic oscillator:

eqn1={x''[t]+ω^2*x[t]==0}

DSolve[eqn1,x[t],t]

The solutions is in the form of a transformation rule. To get the value of x[t] directly use:

x[t]/.DSolve[eqn1,x[t],t]

If there is only one solution, we may also want to get rid of the curly braces by asking for the first component of the solution:

x[t]/.DSolve[eqn1,x[t],t][[1]]

Solve first-order ODE initial value problem for the linear harmonic oscillator:

initial1={x[0]==y0,x'[0]==v0}

sol=x[t]/.DSolve[Join[eqn1,initial1],x[t],t][[1]]

Plot[sol/.{y0->1,v0->2,ω->3},{t,0,10}]

2. Symbolic Solution for Systems of ODE

({{Symbolic Solution for Systems of ODE}, {({{Define the dependent variables, vars={x[t],y[t]}}, {Define ODEs, eqns={eqn1,eqn2}}, {Define initial (or boundary) conditions, inits={x[a]==α,y[a]==α}}, {Give the general solution, sol=vars/.DSolve[eqns,vars,t]}, {Solve an intial value problem , Appart[expr]}, {Plot  x[t] and y[t], Plot[Evaluate[sol],{t,a,b}]}, {Plot phase trajectory , Plot[sol,{t,a,b}]}})}})

vars={Na[t],Nb[t]};

eqns={Na'[t]==-Na[t]/τa,Nb'[t]==Na[t]/τa-Nb[t]/τb}

DSolve[eqns,vars,t]

Ask directly for the values of Na[t] and Nb[t] :

vars/.DSolve[eqns,vars,t][[1]]//Simplify

inits={Na[0]==100,Nb[0]==0};

sol=vars/.DSolve[Join[eqns,inits],vars,t][[1]]

Plot[Evaluate[sol/.{τa->1,τb->0.1}],{t,0,2}]

Special solution is obtained when τa =τb:

Limit[sol,τa->τb]

3. Numerical Solution

({{Numerical Solution for  Single ODE}, {({{Solve the problem numerically, sol=y[t]/.NDsolve[{eqn,conds},y[t],{t,a,b}]}, {Plot the solution, Plot[sol,{t,a,b}]}})}})

Solve numerically first-order ODE for the nonlinear pendulum :

eqn2=θ''[t]+(g/L)*Sin[θ[t]]==0

sol=θ[t]/.NDSolve[{eqn2/.{g->9.81,L->1},θ[0]==Pi/2,θ'[0]==0},θ[t],{t,0,20}][[1]]

Plot[sol,{t,0,20}]

Chaos in Lorentz equations (3 coupled nonlinear ODE describing convetion in the atmosphere):

lorenzEqs = {
		x'[t] == α( -x[t]+y[t]),
		y'[t] == -x[t] z[t] + β x[t] - y[t],
		z'[t] == x[t] y[t] - γ z[t] }

lorenzPlot[values_,ic_,tmax_]:= ParametricPlot3D[Evaluate[{x[t],y[t],z[t]}/.NDSolve[Join[lorenzEqs/.values,ic],{x,y,z},{t,0,tmax},MaxSteps->10000][[1]]],{t,0,tmax},PlotPoints->3000]

A trajectory using the customary parameters illustrates the complexity of this motion:

lorenzPlot[{α->10,β->28,γ->(8/3)},{x[0]==0,y[0]==1,z[0]==-1},50]