Progress in Symbolic Differential Equations

Aram Manaselyan
Wolfram Research, Inc.

(RowDefault[, (Spacer1), , (Spacer1), ])  #WolframTechConf

Outline

Singular solutions for nonlinear ODEs

Improved Method of Undetermined Coefficients

Fractional Differential Equations

Documentation and Tutorials

Singular Solutions for Nonlinear ODEs

By default, DSolve returns a general solution depending on arbitrary parameters for a linear or nonlinear ODE.

For some nonlinear ODEs there can also be singular solutions which cannot be obtained by assigning specific values to the arbitrary constants in the general solution.

To illustrate this, consider the following Clairaut ODE:

DSolve[((y'[x])^2)-y'[x] x+y[x]==0,y[x],x]

The following plot shows the general solution as well as the envelope formed by this family:

Plot[Table[x c-(c^2),{c,-10,10}]//Evaluate,{x,-50,50},PlotRange->All]

The singular solution in this case is the "envelope" y=((x^2)/4):

DSolve[y'[x]^2-y'[x]*x+y[x]==0,y[x],x,IncludeSingularSolutions->True]

Singular Solutions for Nonlinear ODEs

As a second example, let’s consider the following logistic equation:

DSolve[y'[x]==y[x] (1-y[x]),y[x],x]

The following plot shows the general solution as well as the envelope formed by this family:

Plot[Table[((E^x)/((E^x)+(E^c))),{c,-10,10}]//Evaluate,{x,-20,20},PlotRange->All]

The singular solutions in this case are the "equilibrium" solutions y = 0 and y = 1:

DSolve[y'[x]==y[x] (1-y[x]),y[x],x,IncludeSingularSolutions->True]

Improved Method of Undetermined Coefficients

The general solution of a linear inhomogeneous ODE consists of two parts: General solution of homogeneous ODE and particular solution of inhomogeneous equation.

Usually DSolve applies the method of variation of parameters to find the particular solution for a linear inhomogeneous ODE.

This solves an inhomogeneous linear second-order ODE:

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

For some cases the outcome of the method of variation of parameters can be very large and, even unusable.

Improved Method of Undetermined Coefficients

The method of undetermined coefficients can be used for solving inhomogeneous linear ODEs with constant coefficients, if the inhomogeneous term is a polynomial function of independent variable.

sol=DSolve[y''[x]+3y[x]==(x^2)+5x+6, y[x],x]

In version 13.2 we have updated significantly the method of undetermined coefficients, to work also with inhomogeneous terms of the form P (x) (e^(a x)) cos(b x)+Q (x) (e^(a x)) sin(b x) and its sub-forms, where P (x) and Q( x) are polynomial functions of x.

This solves an inhomogeneous linear second-order ODE using the improved method of undetermined coefficients:

DSolve[5y'''[x]+2y[x]==x Exp[-3x] Cos[-2x] ,y[x],x]

Solution obtained by version 13.1:

Improved Method of Undetermined Coefficients

Here is another example of higher order inhomogeneous ODE, solved by the method of undetermined coefficients:

DSolve[Derivative[4][y][x]+5y'''[x]+2y[x]==x (E^(3x)) Cos[3x],y[x],x]

The leaf count of the solution obtained in version 13.1 is 2280.

Improved Method of Undetermined Coefficients

The improved method of undetermined coefficients can also be used for solving inhomogeneous linear systems of ODEs:

(m={{2,-1},{1,-2}})//MatrixForm

g={(2t-5)Cos[3t],(3t+2)Sin[3t]} (E^(-2t));

DSolve[x'[t]-m.x[t]==g,Element[x[t],Vectors[2]],t]

LeafCount[%]

The solution obtained in version 13.1 contains integrals and the leaf count is 852.

Fractional Differential Equations

Fractional differential equations (FDEs) are differential equations involving fractional derivatives (d^α)/d (x^α), where α can be any real positive number.

FDEs are generalizations of the ordinary differential equations.

Currently DSolve can solve only Caputo type FDEs.

Any Caputo-type linear fractional differential equation with constant coefficients can be solved in Wolfram Language.

Example of a Caputo-type FDE:

eqn=CaputoD[y[t],{t,(9/10)}]+10 y[t]==0;

We can solve this FDE using the Laplace transform method:

lteqn=LaplaceTransform[eqn,t,s]/. {y[0]->C[1]}

ltsol=Solve[lteqn,LaplaceTransform[y[t],t,s]]

And finally we apply InverseLaplaceTransform to find the solution of FDE:

InverseLaplaceTransform[ltsol,s,t]

Find the solution directly using DSolve:

sol=DSolve[eqn,y[t],t]

This will verify the solution:

eqn/.sol[[1]]//FullSimplify

Fractional Differential Equations with Initial conditions

For Caputo type FDEs the initial conditions are to be given for integer order derivatives. The number of initial conditions depends on the order α of FDE and is equal to (Ceiling[α]).

This example solves an FDE of order π including four initial conditions:

sol=DSolve[{CaputoD[y[x],{x,Pi}]== y[x],y[0]==1,y'[0]==2,y''[0]==3,y'''[0]==4},y[x],x]

Here is the plot of the solution:

Plot[y[x]/.sol[[1]],{x,0,1}]

Systems of Fractional Differential Equations

DSolve can also handle some linear systems of FDEs with constant coefficients. All equations in the system must have the same order α<1:

eqn1=CaputoD[x[t],{t,(16/17)}]==y[t];
eqn2=CaputoD[y[t],{t,(16/17)}]==-x[t];
ic={x[0]==-3,y[0]==5};

sol=DSolve[Flatten[{eqn1,eqn2,ic}],{x[t],y[t]},t]

Plot the solutions:

Plot[Evaluate[{x[t],y[t]}/.sol[[1]]],{t,0,10}]

Parametrically plot the solution:

ParametricPlot[Evaluate[{x[t],y[t]}/.sol[[1]]],{t,0,20}]

Asymptotic solutions of Fractional Differential Equations

In version 13.2 we have implemented also series solutions of FDEs using AsymptoticDSolveValue.

This finds the asymptotic series solution for a fractional ODE of order 0.7:

asol[n_]:=AsymptoticDSolveValue[{CaputoD[y[x],{x,0.7}]+3 y[x]==0,y[0]==1},y[x],{x,0,n}]

asol[5]

Exact solution can be obtained using DSolveValue:

dsol=DSolveValue[{CaputoD[y[x],{x,0.7}]+3 y[x]==0,y[0]==1},y[x],x]

This compares the asymptotic and exact solutions:

Plot[Evaluate[{dsol,asol[5],asol[10],asol[15],asol[20]}],{x,0,2},PlotLegends->{"Exact","n=5","n=10","n=15","n=20"}]

AsymptoticDSolveValue can be used also for solving FDEs with non constant coefficients:

sol=AsymptoticDSolveValue[{CaputoD[y[x],{x,3/2}]-x y'[x]-y[x]==0,y[0]==1,y'[0]==0},y,{x,0,15}]

Applications of Fractional Differential Equations

Fractional Harmonic Oscillator

This solves the equation of a fractional harmonic oscillator of order 1.9:

sol=DSolveValue[{CaputoD[x[t],{t,1.9}]+ x[t]==0,x[0]==1,x'[0]==0},x[t],t]

Plot the solution:

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

A fractional harmonic oscillator behaves like a damped harmonic oscillator:

dampedOsc=DSolveValue[{x''[t]+a x'[t]+ x[t]==0,x[0]==1,x'[0]==0},x[t],t]//FullSimplify

Fit the data for the fractional oscillator with the solution for the damped oscillator to determine the corresponding value for damping factor a:

data=Cases[InputForm[fracsol],Line[m_]:>m,-1][[1]];
fit=FindFit[data,dampedOsc,a,t]

Compare the results:

Plot[Evaluate[{sol,dampedOsc/.fit}],{t,0,20},PlotLegends->{"Fractional Oscillator","Damped Oscillator"}]

Fractional LC Circuit

This equation describes a fractional LC electrical circuit:

LCeqn[α_,l_,c_]={CaputoD[i[t],{t,α}]+(1/(l c))i[t]==0,i[0]==0.1,i'[0]==0};
LCsol[α_,l_,c_]:=DSolveValue[LCeqn[α,l,c],i[t],t]

We could get the classical solution specifying α=2:

LCsol[2,l,c]

For α=1.9 we have the fractional generalized model for which the solution is:

LCsol[1.9,l,c]

We can plot the solution for various values of equation order α:

Plot[Evaluate[{LCsol[2,1,1],LCsol[1.9,1,1],LCsol[1.8,1,1]}],{t,0,10},PlotLegends->{"α=2","α=1.9","α=1.8"}]

A Dynamical System with Fractional Damping

A. A. Kilbas, H. M. Srivastava and J. J. Trujillo, Theory and Applications of Fractional Differential Equations (Elsevier 2006).

Consider the following dynamical system: a thin plate of mass m and area S that is immersed in a Newtonian fluid of infinite extension with density ρ and viscoelastic constant μ and is connected to a fixed point via a spring with spring constant k. Suppose that the surface of the plate is sufficiently large so that the fluid slows the movement of the plate.

This system can be modeled by a second-order differential equation with fractional damping, using parameters a=((2 S (Sqrt[μ ρ]))/m) and  b=(k/m):

eqn=y''[t]+ a CaputoD[y[t],{t,3/2}]+ b y[t]==f[t];

Solve the model equation for step-like force f(t)=({{{, ({{10, t≤1}, {0, t>1}})}}) with the parameter values a=3, b=1 and plot the solution:

f[t_]=If[t<=1,10,0];
sol=DSolve[{eqn/.{a->3,b->1},y[0]==0,y'[0]==0},y[t],t];
Plot[y[t]/.sol[[1]],{t,0,30}]

Documentation and Tutorials

We have updated the documentation pages for DSolve, DSolveValue and AsymptoticDSolveValue to cover all new features and improvements.

We have updated the tutorial “Symbolic Differential Equation Solving” and added new sections about fractional ODEs, singular solutions etc.

We have updated the tutorial “Working with DSolve: A User's Guide” and added detailed information about recent DSolve features such as Assumptions, IncludeSingularSolutions, Method->Holonomic etc.

Thank you!