Scientific & Engineering Programming

Lecture III

Differential equations with Mathematica

Robert Muszyński

Department of Cybernetics and Robotics
Wrocław University of Science and Technology

Copyright 2017 RM

Derivative notation

Use functions to differentiate expressions or define them

D[]

((dsin x)/(dx)), (((d^2)sin x)/(d(x^2)))

D[Sin[x],x]
D[%,x]

((∂((x^3)+x+(y^4)+(y^2)+(x^3)(y^2)))/(∂x)), (((∂^2)((x^3)+x+(y^4)+(y^2)+(x^3)(y^2)))/(∂(x^2))), (((∂^2)((x^3)+x+(y^4)+(y^2)+(x^3)(y^2)))/(∂x∂y))

D[x^3+x+y^4+y^2+x^3 y^2,x]
D[x^3+x+y^4+y^2+x^3 y^2,{x,2}]
D[x^3+x+y^4+y^2+x^3 y^2,x,y]

((∂f(x,y))/(∂x)), (((∂^2)f(x,y))/(∂(x^2))), (((∂^2)f(x,y))/(∂x∂y))

f[x_,y_]:=x^3+x+y^4+y^2+x^3 y^2;
f[x,y]
D[f[x,y],x]
D[f[x,y],{x,2}]
D[f[x,y],x,y]

d((x^2)+x (y^2)), ((d((x^2)+x (y^2)))/(dx))

Dt[x^2+x y^2,x]    (* total derivative *)
D[x^2+x y^2,x]

Derivative[]

(((∂^2)f(x,y))/(∂(x^2)))=x+(y^2)

Derivative[0,2][u][x_,y_]:=x+y^2

u[x,y]
?u

D[u[a,b],{b,2}]

D[fun[x],x]
%//FullForm

D[fun[x,y],x]
%//FullForm

Apostrophe sign

f'(x)=((df(x))/(dx)), f''(x)=(((d^2)f(x))/(d(x^2)))

f[x_]:=Sin[3x]+x^3
f'[x]         (* Equivalent to D[f[x],x] *)
f''[x]
f''[z]

But be careful!

f'[x^2]   (* compare D[f[x^2],x] *)
Sin'[x]
Sin'[3 x]
Sin'[x^2+3x]
x^3'

g[x_]:=Sin[3x]y+y^3
g[z]
g'[z]

h[x_,y_]:=Sin[3x]y^2+y^3
h[z,v]
h'[z,v]

D[h[z,v],z]
D[h[z,v],z,v]
D[h[z,v],{z,4}]
D[h[z,v],{z,4},{v,2}]

Differentiation

Functions with a single argument

((d(x^n))/(dx)),((df(x))/(dx)),(((d^2)f(x))/(d(x^2)))

D[x^n,x]
D[fun[x],x]
D[fun[x],{x,2}]

((d((Subscript[f, 1])(x)+(Subscript[f, 2])(x)))/(dx)), ((d(Subscript[f, 1])(x)(Subscript[f, 2])(x))/(dx))

D[fun1[x]+fun2[x],x]
D[fun1[x]fun2[x],x]

((df((x^2)))/(dx)), ((d(f^2)(x))/(dx)), ((dx f((x^2)))/(dx))

D[fun[x^2],x]
D[fun[x]^2,x]
D[x fun[x^2],x]
D[%,x]

({{D[f,x], partial derivative (∂/(∂x))f}, {D[f,x,y,…], multiple derivative (∂/(∂x))(∂/(∂y))…f }, {D[f,{x,n}], n(^th) derivative ((∂^n)/(∂(x^n)))f}, {D[f,x,NonConstants->{(Subscript[v, 1]),(Subscript[v, 2]),…}], (∂/(∂x))f with the (Subscript[v, i]) taken to depend on x }})

Functions with several arguments

((∂f(x,y))/(∂x)), ((∂f(x,y))/(∂y)), ((∂f(y,x))/(∂x)), ((df(y(x),x))/(dx)),

expr[x_,y_]:=11y^3+21y^2x+(x+4)^3//Expand
expr[x,y]
D[expr[x,y],x]
D[expr[x,y],y]
D[expr[y,x],x]
D[expr[y[x],x],x]

D[x^2+y^2,{{x,y}}]               (*gradient of the function *)
D[x^2+y^2,{{x,y},2}]         (* Hessian *)
D[{x^2+y^2,x y},{{x,y}}] (* Jacobian *)

Differential equations

Defining

Clear[f];
D[f[t],t]
f'[t]
f'[t]//FullForm

f'(t)+f(t)=1

ode=D[f[t],t]+f[t]==1
ode//TraditionalForm

Solving eqns and utilizing the results

sol=DSolve[ode,f[t],t]

f[t]/.sol[[1]]
f[1]/.sol[[1]]    (* sol contains the rule for f[t] only *)
f[t]/.sol[[1]]/.t->1

Plot[Evaluate[f[t]/.sol[[1]]/.C[1]->1],{t,-5,10}] (* plots work *)

f'[t]/.sol[[1]]  (* but this doesnot work as well *)
f'[1]/.sol[[1]]

g[t_]=f[t]/.sol[[1]];  (* overcome - g[t] evalautes like any other function *)
g[t]
g[1]
g'[t]

Plot[Evaluate[g[t]/.C[1]->1],{t,-5,10}]
Plot[Table[Evaluate[g[t]/.C[1]->param],{param,-5,5}],{t,-5,10}]

f'(t)+f(t)=1, f(0)=2 - initial condition

ode
odei={ode,f[0]==2}  (* condition not necessarily at 0 *)

sol
soli=DSolve[odei,f[t],t]
%//Expand
Plot[Evaluate[f[t]/.%],{t,-5,10}]

y''(x)=y(x), y(0)=1, y'(0)=1 - initial conditions

sol=DSolve[{y''[x]==y[x]},y[x],x]
DSolve[{y''[x]==y[x], y[0]==1},y[x],x]
DSolve[{y''[x]==y[x], y'[0]==1},y[x],x]
DSolve[{y''[x]==y[x], y[0]==1,y'[0]==1},y[x],x]

y[x]/.sol/.{C[1]->1,C[2]->1}

g[x_]=y[x]/.sol[[1]]
Plot[Evaluate[g[x]/.{C[1]->1,C[2]->1}],{x,-10,10}]
gm[x_]=Table[g[x]/.{C[1]->param,C[2]->1},{param,-2,2}]
Plot[Evaluate[gm[x]],{x,-10,10}]

gmm[x_]=Table[g[x]/.{C[1]->param,C[2]->param2},{param,-2,2},{param2,-2,2}]
Plot[Evaluate[gmm[x]],{x,-10,10}]        (* Be careful! 25 functions! *)
g3d[x_,y_]=y[x]/.sol[[1]]/.{C[1]->y,C[2]->y};
Plot3D[Evaluate[g3d[x,y]],{x,-10,10},{y,-10,10},PlotRange->All]

Named versus pure functions

y'(x)+5y(x)=1

ode=y'[x]+5y[x]==1;   (* this we have had above *)
sol = DSolve[ode, y[x],x]
m =sol[[1]]

y[x]/.m
y[1]/.m     (* like previously - no rules for these *)
y'[x]/.m

sol = DSolve[ode, y,x]   (* better approach! *)
m= sol[[1]]     (* the solution in the form of a pure function *)

y[x]/.m
y[1]/.m     (* now it works! *)
y'[x]/.m

ode/.sol   (* to verify the result *)
%//Simplify

Appendix: Pure functions

Clear[x,y,n,m]
Function[x,x^2+4x]
%[n]

Function[{x,y},x^2+4y]
%[n,m]

#^2+4#&
%[n]

#1^2+4#2&
%[n,m]

f[x_]=x^2;                               (* named function *)
Map[f,{2,3,4,5,6}]

Map[#^2&,{2,3,4,5,6}]  (* pure function *)

Multiple solutions

(((y^2)))'(x)=x+11

sol=DSolve[y'[x]^2 == x+11, y,x]

y[x]/.sol[[1]]
y[x]/.sol[[2]]
y[x]/.sol

Plot[Evaluate[y[x]/.sol/.C[1]->10],{x,-15,15}]

Equations with multiple variables

x'(t)=y(t)+z(t)
y'(t)+z(t)-x(t)=0
z'(t}+y(t)=x(t)

sol=DSolve[{x'[t]==y[t]+z[t], y'[t]+z[t]-x[t]==0, z'[t] +y[t]==x[t]}, {x,y,z},t]

{x[t], y[t],z[t]}/.sol[[1]]
traj=%//Simplify  (* but be careful! It may take a lot of time. *)

ParametricPlot3D[Evaluate[traj/.{C[1]->1,C[2]->2,C[3]->3}],{t,0,1}]

Numerical methods for ODEs

y'(x)=y(x),  y(0)=1

DSolve[{y'[x]==y[x],y[0]==1},y,x]

y'(x)=y(x),  y(0)=1
y'(x)=y(x)cos(x+y(x))

(ode={y'[x]==y[x],y[0]==1})//TraditionalForm
DSolve[ode,y,x]
(ode={y'[x]==y[x]Cos[x+y[x]],y[0]==1})//TraditionalForm
DSolve[ode,y,x]

sol=NDSolve[ode,y,{x,0,30}]

y[x]/.sol[[1]]
{y[10],y'[12]}/.sol

Plot[Evaluate[y[x]/.sol],{x,0,30},PlotRange->All]
Plot[Evaluate[{y[x],y'[x],y''[x]}/.sol],{x,0,30},PlotRange->All]

ParametricPlot[Evaluate[{y[x],y'[x]}/.sol],{x,0,30}]

More complicated example - the distribution of heat

Diffusion equation:
((∂w)/(∂t))=α(((∂^2)w)/(∂(x^2)))      ->  w(t,x)

Boundary conditions:
w(0,x)=0; w(t,0)=sin(t); w(t,5)=0

(ode={D[w[t,x],t]==D[w[t,x],x,x],w[0,x]==0,w[t,0]==Sin[t],w[t,5]==0})//TraditionalForm
hsol=NDSolve[ode,w,{t,0,10},{x,0,5}]

Plot3D[Evaluate[w[t,x]/.hsol],{t,0,10},{x,0,5}, PlotRange->All,AxesLabel->{"t","x","T"}]
Plot[Evaluate[w[Pi/2,x]/.hsol],{x,0,5},AxesLabel->{"x","T"}]

Animate[Plot[Evaluate[w[t,x]/.hsol],{x,0,5},AxesLabel->{"x","T"},PlotRange->{-1,1}],{t,0,10},AnimationRunning->False]

Working example: Double pendulum

The system equations

The pendulum state can be described by a vector: ((((Subscript[θ, 1]),(Subscript[θ, 2])))^T) - dependent on time!

The unknown forces: ((((Subscript[F, 1]),(Subscript[F, 2])))^T) - also dependent on time!

Links end points

p1[t_]:={(Subscript[l, 1])Sin[θ1[t]],-(Subscript[l, 1])Cos[θ1[t]]}
p2[t_]:={(Subscript[l, 1])Sin[θ1[t]]+(Subscript[l, 2])Sin[θ2[t]],-(Subscript[l, 1])Cos[θ1[t]]-(Subscript[l, 2])Cos[θ2[t]]}
p1[t]//MatrixForm
p2[t]//MatrixForm

System forces: gravity and tension

vg={0,-g};
fg1=(Subscript[m, 1])vg; fg2=(Subscript[m, 2])vg;
ft1[t_]:=-F1[t]p1[t]/(Subscript[l, 1])+F2[t](p2[t]-p1[t])/(Subscript[l, 2])
ft2[t_]:=-F2[t](p2[t]-p1[t])/(Subscript[l, 2])

By Newton's second law of motion

m p''(t)=∑F

pdeqns={(Subscript[m, 1])p1''[t]==ft1[t]+fg1,   (* pendulum differ. eqns *)
(Subscript[m, 2])p2''[t]==ft2[t]+fg2}

Four equations with variables  (Subscript[θ, 1]),(Subscript[θ, 2]),(Subscript[F, 1]),(Subscript[F, 2]), and their derivatives (but of angles only - thus some equations are algebraic)

The pendulum parameters

pparams={g->9.81,(Subscript[m, 1])->1,(Subscript[m, 2])->1,(Subscript[l, 1])-> 1,(Subscript[l, 2])->1};

The system initial state

pics={θ1[0]==Pi/2,θ2[0]==0,θ1'[0]==0,θ2'[0]==0}

The solution with visualization

psol=NDSolve[{pdeqns,pics}/.pparams,{θ1,θ2,F1,F2},{t,0,15}]

Plot[Evaluate[{θ1[t],θ2[t]}/.psol[[1]]],{t,0,15},PlotLegends->{"(Subscript[θ, 1])","(Subscript[θ, 2])"}]
Plot[Evaluate[{θ1'[t],θ2'[t]}/.psol[[1]]],{t,0,15},PlotLegends->{"(Subscript[θ, 1])'","(Subscript[θ, 2])'"}]

ParametricPlot[Evaluate[p1[t]/.pparams/.psol[[1]]],{t,0,15},AxesLabel->{"x","y"}]
ParametricPlot[Evaluate[p2[t]/.pparams/.psol[[1]]],{t,0,15},AxesLabel->{"x","y"}]

Plot[Evaluate[{F1[t],F2[t]}/.psol[[1]]],{t,0,15},PlotLegends->{"(Subscript[F, 1])","(Subscript[F, 2])"}]

pendul[t_]:={{Line[{{0,0},p1[t],p2[t]}],PointSize[.025],{Red,Point[p1[t]]},{Blue, Point[p2[t]]}}/.psol[[1]]/.pparams}
pendul[0]

Animate[Graphics[pendul[t],PlotRange->{{-2,2},{-2,1}}, Axes->True,Ticks->False,ImageSize->500],{t,0,15,.025},AnimationRate->1/2, AnimationRunning->False]

trace[t_]:={Red,Line[Map[Function[Evaluate[p2[#] /.pparams/.psol[[1]]]],Range[0,t,0.025]]]}
Animate[Graphics[{pendul[t],trace[t]},PlotRange->{{-2,2},{-2,1}}, Axes->True,Ticks->False,ImageSize->500],{t,0,15,.025},AnimationRate->1/2,AnimationRunning->False]

fscale=0.01;
forces[t_]:={{Thick,Red,Arrow[{p1[t],p1[t]-fscale F1[t]p1[t]}],Blue,Arrow[{p2[t],p2[t]-fscale F2[t](p2[t]-p1[t])/(Subscript[l, 2])}]}/. psol[[1]]/.pparams}
forces[1]

Animate[Graphics[{pendul[t],trace[t],forces[t]},PlotRange->{{-2,2},{-2,1}}, Axes->True,Ticks->False,ImageSize->500],{t,0,15,.025},AnimationRate->1/2,AnimationRunning->False]

Alternate approach
(on the basis of Wolfram Double Pendulum example)

The state vector: ((((Subscript[x, 1]),(Subscript[y, 1]),(Subscript[x, 2]),(Subscript[y, 2])))^T) - dependent on time!

The unknown forces: ((((Subscript[λ, 1]),(Subscript[λ, 2])))^T) - dependent on time! ((Subscript[F, 1]) and (Subscript[F, 2]) on the figure and previously)

Again by Newton's second law of motion -  m (x')'(t)==∑(Subscript[λ, x]) , m (y')'(t)==∑(Subscript[λ, y])

npdeqns={(Subscript[m, 1])x1''[t]==λ1[t]((x1[t])/(Subscript[l, 1]))-λ2[t]((x2[t]-x1[t])/(Subscript[l, 2])), (Subscript[m, 1])y1''[t]==λ1[t] ((y1[t])/(Subscript[l, 1]))-λ2[t]((y2[t]-y1[t])/(Subscript[l, 2]))-(Subscript[m, 1]) g,
(Subscript[m, 2])x2''[t]==λ2[t]((x2[t]-x1[t])/(Subscript[l, 2])),
(Subscript[m, 2])y2''[t]==λ2[t] ((y2[t]-y1[t])/(Subscript[l, 2]))- (Subscript[m, 2])g};

The constraints reflecting the constant length of the pendulum rods

npaeqns={((x1[t])^2)+((y1[t])^2)==((Subscript[l, 1])^2),(((x2[t]-x1[t]))^2)+(((y2[t]-y1[t]))^2) ==((Subscript[l, 2])^2)};

The system initial state

pics     (* we have had in the previous case *)

npics={x1[0]==1,y1[0]==0,x2[0]==1,y2[0]==-1,y1'[0]==0,y2'[0]==0}

npsol=NDSolve[{npdeqns,npaeqns,npics}/.pparams,{x1,y1,x2,y2,λ1,λ2},{t,0,15},Method->{"IndexReduction"->{"Pantelides", "ConstraintMethod"->"Projection"}}]//First

Plot[Evaluate[{x1[t],x2[t]}/.npsol],{t,0,15},PlotLegends->{"(Subscript[x, 1])","(Subscript[x, 2])"}]
Plot[Evaluate[{x1'[t],x2'[t]}/.npsol],{t,0,15},PlotLegends->{"(Subscript[x, 1])'","(Subscript[x, 2])'"}]

ParametricPlot[Evaluate[{x1[t],y1[t]}/.npsol],{t,0,15}]
ParametricPlot[Evaluate[{x2[t],y2[t]}/.npsol],{t,0,15}]

Plot[Evaluate[{λ1[t],λ2[t]}/.npsol],{t,0,15},PlotLegends->{"(Subscript[λ, 1])","(Subscript[λ, 2])"}]

?pendul   (* we had the pendulum visualization *)

npendul[t_]:={{Line[{{0,0},{x1[t],y1[t]},{x2[t],y2[t]}}],PointSize[.025],{Red,Point[{x1[t],y1[t]}]},{Blue, Point[{x2[t],y2[t]}]}}/. npsol}
pendul[0]

Animate[Graphics[npendul[t],PlotRange->{{-2,2},{-2,1}}, Axes->True,Ticks->False,ImageSize->500],{t,0,15,.025},AnimationRate->1/2, AnimationRunning->False]

ntrace[t_]:={Blue,Line[Map[Function[Evaluate[{x2[#],y2[#]} /.npsol]],Range[0,t,0.025]]]}
trace[0.1]

Animate[Graphics[{npendul[t],ntrace[t]},PlotRange->{{-2,2},{-2,1}}, Axes->True,Ticks->False,ImageSize->500],{t,0,15,.025},AnimationRate->1/2,AnimationRunning->False]

Comparison

Both models on one plot

Animate[Graphics[{npendul[t],pendul[t],ntrace[t],trace[t]},PlotRange->{{-2,2},{-2,1}}, Axes->True,Ticks->False,ImageSize->500],{t,0,15,.025},AnimationRate->1/2,AnimationRunning->False]

Differences - what are the reasons?

Let’s play with initial conditions: the pendulum hanging down

The original model

pics={θ1[0]==0,θ2[0]==0,θ1'[0]==0,θ2'[0]==0};
psol=NDSolve[{pdeqns,pics}/.pparams,{θ1,θ2,F1,F2},{t,0,15}];
Animate[Graphics[{pendul[t]},PlotRange->{{-2,2},{-2,1}}, Axes->True,Ticks->False,ImageSize->500],{t,0,15,.025},AnimationRate->1/2,AnimationRunning->False]

The alternate model

npics={x1[0]==0,y1[0]==-1,y1'[0]==0,x2[0]==0,y2[0]==-2,y2'[0]==0};
npsol=NDSolve[{npdeqns,npaeqns,npics}/.pparams,{x1,y1,x2,y2,F1,F2},{t,0,15},Method->{"IndexReduction"->{"Pantelides", "ConstraintMethod"->"Projection"}}];
Animate[Graphics[{npendul[t]},PlotRange->{{-2,2},{-2,1}}, Axes->True,Ticks->False,ImageSize->500],{t,0,15,.025},AnimationRate->1/2,AnimationRunning->False]

The state variables

Plot[Evaluate[{θ1[t],θ2[t]}/.psol[[1]]],{t,0,15},PlotLegends->{"(Subscript[θ, 1])","(Subscript[θ, 2])"}]
Plot[Evaluate[{x1[t],x2[t]}/.npsol[[1]]],{t,0,15},PlotLegends->{"(Subscript[x, 1])","(Subscript[x, 2])"}]

The pendulum placed sideways from its resting positon

Original model

tmax=1500;   (* see disturbed (10^-10) for shorter time (15) *)
pics={θ1[0]==Pi,θ2[0]==Pi,θ1'[0]==0,θ2'[0]==0};
psol=NDSolve[{pdeqns,pics}/.pparams,{θ1,θ2,F1,F2},{t,0,tmax}];
Plot[Evaluate[{θ1[t],θ2[t]}/.psol[[1]]],{t,0,tmax},PlotLegends->{"(Subscript[θ, 1])","(Subscript[θ, 2])"}]

tmax=15;
npics={x1[0]==0,y1[0]==1,y1'[0]==0,x2[0]==0,y2[0]==2,y2'[0]==0};
psol=NDSolve[{npdeqns,npaeqns,npics}/.pparams,{x1,y1,x2,y2,F1,F2},{t,0,15},Method->{"IndexReduction"->{"Pantelides", "ConstraintMethod"->"Projection"}}];
Plot[Evaluate[{x1[t],x2[t]}/.psol[[1]]],{t,0,tmax},PlotLegends->{"(Subscript[x, 1])","(Subscript[x, 2])"}]

Rules transforming the original equations to the alternate ones

rules={x1[t]->(Subscript[l, 1])Sin[θ1[t]],y1[t]->-(Subscript[l, 1]) Cos[θ1[t]],x2[t]->(Subscript[l, 1]) Sin[θ1[t]]+(Subscript[l, 2]) Sin[θ2[t]],y2[t]->-(Subscript[l, 1]) Cos[θ1[t]]-(Subscript[l, 2]) Cos[θ2[t]]}

npdeqns/.rules
npdeqns/.rules/.D[rules,{t,2}]

Another example: RLC circuit

The system equations
(see Applications in NDSolve help)

By component laws and Kirchoff’s laws:

components={r iR[t]==vR[t],l iL'[t]==vL[t],iC[t]==c vC'[t]};
connections={vR[t]+vL[t]+vC[t]-v1[t]==0,iL[t]==iR[t],iR[t]==iC[t]};

The parameters and initial conditions

params={r->5,l->10^-2,c->10^-4};ic={iR[0]==0,iL[0]==0,iC[0]==0,vR[0]==0,vL[0]==0,vC[0]==0};

A step response

tmax=0.1;
v1[t_]:=UnitStep[t-.01];
csol=NDSolve[{components,connections,ic}/.params,vC,{t,0,tmax}];
Plot[Evaluate[{v1[t],vC[t]}/.csol],{t,0,tmax}]
Plot[Evaluate[{v1[t],vC[t]}/.csol],{t,0,tmax},PlotRange->All]