What’s New in DSolve

Devendra Kapadia

Wolfram Research, Inc.

Outline of the Talk

In this talk, I will give a fairly complete introduction to the following new DSolve features in Mathematica 10:

Improvements for initial-value problems.

Solutions of delay differential equations.

Differential equations with events.

Using DSolveValue to obtain solution expressions.

Initial Value Problems

The General Algorithm for IVPs

Given an initial value problem, DSolve first attempts to find the general solution for the ODE and then fits the arbitrary constants to the initial conditions.

A Linear Second-Order ODE

Consider the following IVP for a linear second-order ODE:

linearivp={y''[x]-3 y'[x]+2 y[x]==0,y[0]==1,y'[0]==1};

DSolve returns a particular solution in this case:

DSolve[linearivp,y[x],x]

Internally, we first find the general solution:

sol=DSolve[y''[x]-3 y'[x]+2 y[x]==0,y,x]

Next, we apply the initial conditions:

init={(y[0]/. sol[[1]])==1,(y'[0]/. sol[[1]])==1}

We now solve this system:

Solve[init,{C[1],C[2]}]

Finally, the values of arbitrary constants are plugged into the general solution:

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

This scheme has limitations for nonlinear ODEs because the general solution have a complicated expresssion given in terms of special functions.

The Pendulum Equation

The  general solution for the pendulum equation depends on JacobiAmplitude.

pendulumeqn={θ''[t]==-Sin[θ[t]]};

Quiet[sol=DSolve[pendulumeqn,θ,t]]

FullSimplify[pendulumeqn/. sol]

The equation for the simple pendulum does not depend on the first derivative. Hence, we apply a special (two-step) method for solving the corresponding initial value problem.

pendulumivp={θ''[t]==-Sin[θ[t]],θ[0]==0,θ'[0]==1};

sol=DSolve[pendulumivp,θ,t]

Here is a plot of the solution:

Plot[θ[t]/. sol,{t,0,20},Filling->Axis,FillingStyle->LightBlue,PlotStyle->{Thickness[0.006],Red}]

Solutions in terms of Hypergeometric Functions

This initial value problem is solved in terms of Hypergeometric2F1 in Mathematica 10.

eqn={y''[x]==4 y'[x]+(1/((y'[x])^2)),y[1]==3,y'[1]==1};

sol=DSolve[eqn,y,x]

The solution may be verified as follows.

FullSimplify[eqn/. sol]

Delay Differential Equations

A First-Order Delay Equation

Here is a first-order delay equation which is solved using the method of steps on the interval {t, 0,2};

sol=DSolve[{x'[t]==((x[t-1])^2),x[t/;t<0]==a t},x,{t,0,2}]

The solution may be plotted for different values of ‘a’ as follows:

Plot[Table[x[t]/. sol[[1]],{a,1,5}],{t,-1,2},PlotStyle->Red]

A System of Delay Equations

The following solves a system of delay differential equations:

eqns={x'[t]==a y[t-1]+y[t-3],y'[t]==x[t-1],x[t/;t≤0]==t,y[t/;t≤0]==(t^2)};

sol=DSolve[eqns,{x,y},{t,0,5}]

Here is a plot of the solution:

ParametricPlot[Evaluate[Table[{x[t],y[t]}/. sol,{a,-1,2,(1/3)}]],{t,-1,5},WorkingPrecision->20,PlotStyle->Red,Exclusions->None]

Chaotic Behavior

The solution of the following DDE exhibits chaotic behavior:

sol=DSolve[{x'[t]==Sign[x[t-2]]-x[t-2],x[t/;t≤0]==0.1},x,{t,0,20}];

Plot[x[t]/. sol[[1]],{t,0,20},PlotStyle->{Thickness[0.005],Red}]

ODEs with Events

An ODE with a Time-Dependent Event

Here is an ODE with a time-dependent event specified by WhenEvent:

dsol=DSolve[{x'[t]==x[t],x[0]==1,WhenEvent[Mod[t,1]==0,x[t]->x[t]+10]},x,{t,0,3}]

A plot of the solution clearly shows the jumps at t = 1 and t = 2:

Plot[x[t]/. dsol[[1]],{t,0,3},Exclusions->None,PlotStyle->Red]

The Language of Events

Solve an ODE with a state-dependent event.

sol=DSolve[{y''[t]==-(981/100),y[0]==5,y'[0]==0,WhenEvent[y[t]==0,y'[t]->-(95/100) y'[t]]},y,{t,0,10}]

This system models a bouncing ball.

Plot[y[t]/. sol,{t,0,10},Filling->Axis]

Stop the integration of an ODE when an event occurs.

DSolve[{x'[t]==-x[t],x[0]==1,WhenEvent[t==3,"StopIntegration"]},x,{t,0,4}]

Remove an event after it has occurred once.

sol=DSolve[{x''[t]+x[t]==0,x[0]==0,x'[0]==1,WhenEvent[x[t]==0,{x'[t]->-x'[t],"RemoveEvent"}]},x,{t,0,14}]

Plot[x[t]/. sol,{t,0,14}]

Ball Bouncing Down Steps

The following models a ball bouncing down steps using DiscreteVariables to model the steps.

sol=DSolve[{y''[t]==-9.8,y[0]==13.5,y'[0]==5,a[0]==13.5,WhenEvent[y[t]==a[t],y'[t]->-0.75 y'[t]],WhenEvent[Mod[t,1]==0,a[t]->a[t]-1]},{y,a},{t,0,10},DiscreteVariables->a]

Plot[Evaluate[{y[t],a[t]}/. sol],{t,0,10},Filling->{2->0},Exclusions->None,ImageSize->Medium]

Ball in a Square Box

In a square box, model a ball that changes direction upon impact with the side walls:

sol=DSolve[{x'[t]==a[t],y'[t]==b[t],x[0]==0,y[0]==0,a[0]==1,b[0]==Rationalize[(Sqrt[2]),.01],WhenEvent[((x[t])^2)==1,a[t]->-a[t]],WhenEvent[((y[t])^2)==1,b[t]->-b[t]]},{x,y},{t,0,100},DiscreteVariables->{a,b}];

LeafCount[sol]

ParametricPlot[{x[t],y[t]}/. sol,{t,0,100},Frame->True,FrameTicks->None,PlotRange->1,Axes->False,PlotPoints->200]

With an irrational initial velocity vector, there are no longer periodic solutions:

sol=DSolve[{x'[t]==a[t],y'[t]==b[t],x[0]==0,y[0]==0,a[0]==1,b[0]==(Sqrt[2]),WhenEvent[((x[t])^2)-1,a[t]->-a[t]],WhenEvent[((y[t])^2)-1,b[t]->-b[t]]},{x,y},{t,0,100},DiscreteVariables->{a,b}];

LeafCount[sol]

ParametricPlot[{x[t],y[t]}/. sol,{t,0,100},Frame->True,FrameTicks->None,PlotRange->1,Axes->False,PlotPoints->200]

A Proportional-Derivative Controller

Model the position of a moving body with 1 kg mass:

system={x''[t]==u[t],x'[0]==x[0]==0,u[0]==1};

Use a sampled proportional-derivative (PD) controller to keep the position constant:

kp=1;td=1;xref=1;
control=WhenEvent[Mod[t,1]==0,u[t]->kp (xref-x[t]-td x'[t])];

Simulate and visualize:

sol=DSolve[{system,control},{x,u},{t,0,15},DiscreteVariables->u]

Plot[Evaluate[{xref,x[t],u[t]}/. sol],{t,0,15},Exclusions->None]

DSolveValue

Basic Uses

DSolveValue returns an expression for the symbolic solution of a differential equation:

DSolveValue[{y'[x]==2 x,y[0]==0},y[x],x]

A pure function is returned if we use y instead of y[x] in the second argument:

sol=DSolveValue[{y'[x]==2 x,y[0]==0},y,x]

Plot[sol[x],{x,0,10}]

Using Arbitrary Expressions in DSolveValue

The following returns the general solution for the Airy equation:

DSolveValue[{y''[x]==x y[x]},y[x],x]

The derivative of the general solution can be obtained as follows:

DSolveValue[{y''[x]==x y[x]},y'[x],x]

Find the value of the solution for an ODE at a point:

DSolveValue[{y''[x]+y[x]==0,y[0]==0,y'[0]==1},y[(π/3)],x]

Elliptic Wave Rectification

The following models a rectified elliptic wave:

Quiet[sol=DSolveValue[{y''[t]==-Sin[y[t]],y[0]==0,y'[0]==1,WhenEvent[y[t]==0.,{y[t]->0,y'[t]->-y'[t]}]},y,{t,0,12}]]

Plot[sol[t],{t,0,10.2},PlotStyle->Red,Exclusions->None,WorkingPrecision->20,Filling->Axis,FillingStyle->LightBlue]

Conclusion

There are significant improvements for the functionality and interface of DSolve in Mathematica 10.

We hope to continue with other improvements (e. g. systems of ODEs, PDEs) in the near future.