Chapter 3
Differential equations, DSolve related

3.1 How to parse an ODE?

This code takes one ode, using same input as given to Mathematica’s DSolve and verifies it is valid. I use this in my ode solver. This is the very first function called to verify the input is valid.

 
symbolQ[_Symbol] = True; 
symbolQ[_]       = False; 
has[e_,h_]       := Not[FreeQ[e,h]] 
 
(*This function thanks to Carl Woll,see https://mathematica.stackexchange.com/questions/151850/using-cases-and-when-to-make-input-a-list-or-not*) 
getPatterns[expr_,pat_]:=Last@Reap[expr/. a:pat:>Sow[a],_,Sequence@@#2&]; 
 
(*Version Dec 2, 2024 by Nasser M. Abbasi*) 
(*Takes ode and IC, using same syntax as Mathematica's DSolve, and verifies it is valid*) 
(*Returns list of ode,IC,ode order if valid, else it generates Abort message*) 
 
parseOneOde[odeAndIc_,y_Symbol[x_Symbol],x_Symbol]:=Module[{ode,ic,constX,alg,LHS,z,n,indeps,deps,der,odeOrder}, 
 
If[Not[symbolQ[y]], 
    Print["Invalid dependent variable given ",y]; 
    Abort[] 
]; 
 
If[Not[symbolQ[x]], 
     Print["Invalid independent variable given ",x]; 
     Abort[] 
]; 
 
If[Head[odeAndIc]===List, 
    {constX,ic,alg,ode}=Internal`ProcessEquations`SeparateEquations[Flatten[odeAndIc],{x},{y}] 
    , 
    {constX,ic,alg,ode}=Internal`ProcessEquations`SeparateEquations[{odeAndIc},{x},{y}] 
]; 
 
If[Length[Flatten[constX]]!=0, 
     Print["Does not support constraints on independent variable ",constX]; 
     Abort[] 
]; 
 
If[Length[Flatten[alg]]!=0, 
    Print["Does not support algebratic equations ",alg]; 
    Abort[] 
]; 
 
If[Length[ode]>1, 
    Print["Does not support more than one ode at this time",ode]; 
    Abort[] 
]; 
 
If[Length[ode]==0, 
    Print["No ode in ",y[x]," found"]; 
    Abort[] 
]; 
 
(*extract the ode from the list now we know there is only one*) 
ode=First[Flatten@ode]; 
 
(*this makes it easier to parse*) 
LHS=First@ode-Last@ode; 
 
If[FreeQ[LHS,y], 
     Print["Error: ode ",ode," has no ",y]; 
     Abort[] 
]; 
 
If[FreeQ[LHS,x], 
    Print["Error. ode ",ode," has no ",x]; 
    Abort[] 
]; 
 
deps=getPatterns[LHS,y[]]; 
 
If[deps=!={}, 
     Print["makeODE:: Can not have ",y[]," in ",ode]; 
     Abort[] 
]; 
 
deps=getPatterns[LHS,y[z_.]]; 
 
(*this checks there is no y on its own*) 
If[Length[Cases[Variables[LHS],Verbatim[y]]]!=0, 
     Print["makeODE:: Can not have ",y," with no argument ",x," in ode ",ode]; 
     Abort[] 
]; 
 
(*this code checks all the independent variable in y[x] are x*) 
If[deps=!={}, 
     deps=Union[Cases[deps,y[z_]:>z]];(*this extracts all the independent variables*) 
     If[Length[deps]>1, 
         Print["makeODE:: independent variable must be ",x," in ",ode]; 
         Abort[] 
     ]; 
 
     If[First[deps]=!=x, 
         Print["makeODE:: independent variable must be ",x," in ",ode]; 
         Abort[] 
     ] 
]; 
 
der = getPatterns[LHS,Derivative[n_][y][x]]; 
 
If[der==={}, 
    Print["makeODE::No differential equation found in ",ode]; 
    Abort[] 
]; 
 
der = getPatterns[LHS,Derivative[n_][y][_]]; 
 
indeps = Union@Cases[der,Derivative[n_][y][z_]:>z]; 
 
If[Or[Length@indeps>1,First@indeps=!=x], 
     Print["Found wrong independent variable in derivative in ",ode]; 
     Abort[] 
]; 
 
ic = Union@ic; 
 
odeOrder = First@Flatten@Internal`ProcessEquations`DifferentialOrder[ode,{x},{y}]; 
 
If[Length@ic>0, 
    If[Not[AllTrue[ic,Head[#]===Equal&]], 
         Print["Invalid initial conditions found in ",ic]; 
         Abort[] 
    ]; 
 
    If[Length@ic>odeOrder, 
         Print["Too many initial conditions given for the order of the ode ",ic]; 
         Abort[] 
    ] 
]; 
 
Print["Successful Parsing"]; 
Print["ode=",ode]; 
Print["y=",y]; 
Print["x=",x]; 
Print["ode order=",odeOrder]; 
Print["IC=",ic]; 
 
{ode,ic,odeOrder} 
]
 

These are examples using the above function

parseOneOde[{y''[x] == y[x]*Sin[x], y[0] == 0, y'[0] == 1}, y[x], x] 
 
    {y''[x] == Sin[x] y[x], {y[0] == 0, y''[0] == 1}, 2}
 

parseOneOde[{y''[x] == y[x]*Sin[x], y[x] == 0, y'[0] == 1}, y[x], x] 
 
    "Does not support algebratic equations ", {y[x] == 0} 
    $Aborted
 

parseOneOde[{y''[x] == y[]*Sin[x], y[Pi] == 0, y'[0] == 1}, y[x], x] 
 
    makeODE:: Can not have y[] in y''[x]==Sin[x] y[] 
    $Aborted
 

parseOneOde[{Sin[x] y''[x] + Derivative[y[x], {x, 7}] == y[x]*Sin[x] - y[x], y[Pi] == 0, y'[0] == 1}, y[x], x] 
 
   {Derivative[y[x],{x,7}]+Sin[x] y''[x]==-y[x]+Sin[x] y[x],{y[Pi]==0,y'[0]==1},2}
 

parseOneOde[{Sin[x] y''[x] == y[z]*Sin[x] - y[x], y'[0] == 1}, y[x],  x] 
 
   makeODE:: independent variable must be , x, in , Sin[x] Derivative[2][y][x] == -y[x] + Sin[x] y[z]] 
   $Aborted
 

3.2 How to find the dependent variable of an ode?

Given a single ode, how to find the dependent variable?

Internal`ProcessEquations`FindDependentVariables[y''[x] + y[x]*y'[x] == x, x] 
 
   {y}
 

Internal`ProcessEquations`FindDependentVariables[x'[t] == 0, t] 
 
   {x}
 

Notice the second argument must be the independent variable.

3.3 How to find if ode is linear or not?

There is no builtin function I could find in Mathematica to do this. But one way is the following. We will have a function that returns the coefficients of the ode and the RHS. i.e. given an ode \(A y''+ By' + C y=f(x)\), it will return a list \(\{A,B,C,f\}\). Next to check if the ode is linear or not, all what we have to do is check that \(A,B,C\) do not have \(y\) in them. This include \(y,y',y''\) and so on.

So for a first order ode, the result of calling this function, which we will get getODEcoeff is always a list. The Length of the list is 2 more than the order of the ode. If a coefficient is missing, zero is placed in that slot. For example, for 3 y''[x] + x y[x]==Sin[x] then getODEcoeff[ode,y[x],x] will return {3,0,x,Sin[x] and for x y''''[x] + x y[x]==Sin[x] it will return {x,0,0,0,x,Sin[x]} and for first order ode y'[x]+f[x]*y[x]==9 it will return {1,f[x],9}.

Now we simply check that all element of this list have no \(y\) in them. Here is the function

(to do)

3.4 DSolveIntegrals package

Where did I get this from?

"Mathematica can handle partial differential equations via the DSolveIntegrals package. These arise in chemical contexts in the 1D wave equation, 3D wave equation, 3D diffusion equation, Time-dependent and Time independent Schrödinger equation. Hermite showed that the quintic equation could be solved by elliptic functions"

3.5 Extracting DSolve solutions

one way

Clear[y, x] 
DSolve[{y'[x]^2 == 1 - y[x]^2, y[0] == 0}, y[x], x] 
 
{{y[x] -> -Sin[x]}, {y[x] -> Sin[x]}}
 

And now

\verb|y[x] /. %| 
 
   {-Sin[x],Sin[x]}
 

3.6 How to solve ODE using power series method?

For an ode, that has just ordinary point at \(x=x_0\), we can find power series solution near \(x_0\) as follows. Assume the ode is

\[ u''(t)+\frac {1}{10} u(t)^3 u'(t)^2+4 u(t)=0 \]

with intitial conditions \(u(0)=1,u'(0)=1\), then

findSeriesSolution[t_, nTerms_] := Module[{pt = 0, u, ode, s0, s1, ic, eq, sol}, 
  ic = {u[0] -> 1, u'[0] -> 1}; 
  ode = u''[t] + 4 u[t] + 1/10 u[t]^3 u'[t]^2; 
  s0 = Series[ode, {t, pt, nTerms}]; 
  s0 = s0 /. ic; 
  roots = Solve@LogicalExpand[s0 == 0]; 
  s1 = Series[u[t], {t, pt, nTerms + 2}]; 
  sol = Normal[s1] /. ic /. roots[[1]] 
  ]
 

and now call it with

seriesSol = findSeriesSolution[t, 6]
 

It returns

\[ \frac {445409479 t^8}{840000000}+\frac {8878343 t^7}{10500000}-\frac {277427 t^6}{600000}-\frac {12569 t^5}{50000}+\frac {1607 t^4}{2000}-\frac {29 t^3}{50}-\frac {41 t^2}{20}+t+1 \]

Update: As of recent versions, we can just do this

ode = u''[t] + 4 u[t] + 1/10 u[t]^3*u'[t]^2 == 0; 
ic = {u[0] == 1, u'[0] == 1}; 
AsymptoticDSolveValue[{ode, ic}, u[t], {t, 0, 6}]
 

\[ -\frac {277427 t^6}{600000}-\frac {12569 t^5}{50000}+\frac {1607 t^4}{2000}-\frac {29 t^3}{50}-\frac {41 t^2}{20}+t+1 \]

3.7 Replacing \(y,y',y''\) in an equation

Suppose we have \(u''(t)+u'(t)+u(t)=3 \cos (2 t)\) and we wanted to find a particula solution by replacing \(u\) in the differential equation by some guess for a particular solution. Then do

ode = Derivative[2][u][t] + Derivative[1][u][t] + u[t] == 3*Cos[2*t]; 
ode /. u -> (c1*Cos[#1] + c2*Sin[#1] & ) 
 
Out[2]= c2*Cos[t] - c1*Sin[t] == 3*Cos[2*t]
 

3.8 How to classify singular points for ODE?

Given an ode, such as \(x(1-x) y''(x)+(c-(a+b+1)x)y'(x)-a b y(x)=0\) and we want to classify the singular points (finite and infinite).

We write it as \(y''(x)+p(x) y'(x)+q(x)=0\) and then follow standard procedures. In this example, we want to classify points \(0,1,\infty \).

This small function I wrote for a HW will do this. Pass it \(p(x),q(x)\) and list of the points to classify.

checkForSingularity[p_, q_, at_List, x_] := Module[{p0, q0, t, r, r0}, 
  r = First@Last@Reap@Do[ 
       If[at[[n]] === Infinity, 
        p0 = p /. x -> (1/t); 
        p0 = (2 t - p0)/t^2; 
        q0 = (q /. x -> (1/t))/t^4; 
        r0 = checkForSingularity[p0, q0, {0}, t]; 
        Sow[Flatten[{Infinity, Rest[Flatten@r0]}]]; 
        , 
        r0 = {at[[n]], Limit[(x - at[[n]])  p, x -> at[[n]]], 
          Limit[(x - at[[n]])^2 q, x -> at[[n]]]}; 
        Sow@r0 
        ] 
       , {n, 1, Length@at} 
       ]; 
  r 
  ]
 

To use it on the above example

ClearAll[c, a, b, x]; 
m = checkForSingularity[(c - (a + b + 1) x)/(x (1 - 
        x)), (-a b)/(x (1 - x)), {0, 1, Infinity}, x]; 
Grid[Join[{{"point", "limit x p(x)", "limit x^2 q(x)"}}, m], 
 Frame -> All]
 

Gives

{{"point", "limit x p(x)", "limit x^2 q(x)"}, 
 {0,            c,                 0}, 
 {1,        1 + a + b - c,         0}, 
 {\[Infinity], 1 - a - b,        a b} 
}
 

Since all points have finite limits, then all \(0,1,\infty \) are removable singularities.

I did not test this function too much, but it works for the HW I did :)

3.9 How to replace \(y(x)\) by \(e^{s(x)}\) in an ode?

Given ode \(y^{'''}=x y(x)\) we want to replace \(y(x)\) by \(e^{s(x)}\)

ClearAll[y,x,s] 
eq = y'''[x] == x y[x] 
eq = eq /. y -> Function[{x}, Exp[s[x]]] 
Simplify@Thread[eq/Exp[s[x]], Equal]
 

Gives \(x=s^{'''}(x)+s'(x)^3+3 s'(x) s''(x)\)

3.10 rewrite ODE as \(y''=\text {RHS}\)

Given y'[x]+3 Sin[x]-3-Cos[x]-2 y[x]==0 how to rewrite it to be y'[x]==3+Cos[x]-3 Sin[x]+2 y[x] i.e. put y'[x] on one side, and the rest on the other side?

ClearAll[f,x,y] 
expr=y'[x]+3 Sin[x]-3-Cos[x]-2 y[x]==0; 
lhs=expr/.lhs_==rhs_:>lhs; 
rhs=expr/.lhs_==rhs_:>rhs; 
rhs=-(lhs/.(y'[x]+any_.):>any)+rhs; 
expr=y'[x]==rhs
 

This makes it easier to make it in form \(y'(x)=f(x,y)\).

3.11 How to trace DSolve running?

Use this

Block[{DSolve`print=Print}, 
    DSolve[{ode,ic},y[t],t] 
]
 

3.12 How to put time out on integrate in DSolve?

see https://mathematica.stackexchange.com/questions/120364/why-cant-dsolve-find-a-solution-for-this-ode/120650#120650 by Michael E2.

ClearAll[withTimedIntegrate]; 
SetAttributes[withTimedIntegrate, HoldFirst]; 
withTimedIntegrate[code_, tc_] := Module[{$in}, 
   Internal`InheritedBlock[{Integrate}, 
    Unprotect[Integrate]; 
    i : Integrate[___] /; ! TrueQ[$in] := 
     Block[{$in = True}, 
      TimeConstrained[i, tc, Inactivate[i, Integrate]] 
      ]; 
    Protect[Integrate]; 
    code 
    ] 
   ]; 
 
withTimedIntegrate[{dsol} = DSolve[ode == 0, y, x], 1]; // AbsoluteTiming 
dsol
 

3.13 How to find the order of an ODE?

Use

Clear["Global`*"]; 
ode=y''[x]+3*y'[x]+3==0; 
Internal`ProcessEquations`DifferentialOrder[ode,{x},{y}] 
 
(*   {{2}} *)
 

So it is order 2. Notice the input must be lists [ode,{x},{y}], it will hang if you use [ode,x,y]

3.14 How to find the degree of an ODE?

Degree of an ode is the degree of the highest derivative in the oder. This needs a helper function

(*This function thanks to  Carl Woll, see https://mathematica.stackexchange.com/questions/151850/using-cases-and-when-to-make-input-a-list-or-not*) 
Clear["Global`*"]; 
getPatterns[expr_, pat_] :=   Last@Reap[expr /. a : pat :> Sow[a], _, Sequence @@ #2 &];
 

And now do

getDegreeOfOde[ode_,y_,x_]:=Module[{maxDer,p,d,der}, 
    der    = getPatterns[ode,Power[Derivative[_.][y][x],_.]] ; 
    p      = Flatten[Internal`ProcessEquations`DifferentialOrder[#,{x},{y}]&/@der]; 
    maxDer = Extract[der,Position[p,Max[p]]]; 
    Abs[First[maxDer/.(Derivative[_.][y][x])^(d_.):>d]] 
]
 

Examples of usage

ode=y'''[x]+y''[x]^4+3*y'[x]^8+3*y[x]^8==0; 
getDegreeOfOde[ode,y,x] 
   (*1*) 
 
ode=y''[x]^4+3*y'[x]^8+3*y[x]^8==0; 
getDegreeOfOde[ode,y,x] 
   (*4*) 
 
ode=y''[x]+3*y'[x]^8+3*y[x]^8==0; 
getDegreeOfOde[ode,y,x] 
   (*1*)
 

3.15 How to move all terms with y to one side of ode?

Given an ode (or any equation), and we want to move all terms with \(y(x)\) to left side and everything else to right side. This makes it easier to see the forcing function. Select is used for this

moveAllYToOneSide[ode_Equal,y_Symbol,x_Symbol]:=Module[{expr=ode,termsWithNoY,termsWithY}, 
expr=expr[[1]]-expr[[2]]; 
termsWithNoY=Select[expr,FreeQ[#,y]&]; 
termsWithY=Select[expr,Not[FreeQ[#,y]]&]; 
expr=termsWithY==-termsWithNoY 
]
 

Call it as

ode=y[x]*y'[x]+Sin[x]+3-1/y[x]==Sin[y[x]]+Pi*x*y[x]; 
moveAllYToOneSide[ode ,y,x]
 

Gives

-Sin[y[x]]-1/y[x]-Pi x y[x]+y[x] y'[x]==-3-Sin[x]