Notes

This notebook contains ‘module’ that allows to integrate ODE along prescried paths in complex domain. The paths are drawed on the same plot so it is possible to see the sketch of the graph of solution.

Author: Radosław Antoni Kycia   ( rkycia@mimuw.edu.pl )

Support: This project was supported by the Warsaw Center of Mathematics and Computer Science from the funds of the Polish Leading National Research Centre (KNOW).

Legal Notica: The  notebook/library is in ‘under developement’ stage. It should work, however if it is not a case please e-mail Author. Author is not responisble for any dameges, data losses which results form the use of this notebook. Use it at your own risk. 

The notebook is distributed under GNU GENERAL PUBLIC LICENSE version 3:
http://www.gnu.org/copyleft/gpl.html

Request: If you find this notebook somehow usefull then you are kindly asked to cite paper which describes these ideas:

Radosław A. Kycia ‘Scheme for numerical investigation of movable singularities of the complex valued solutions of ordinary differential equations’, accepted. 

The references therein can be also useful for You to develop your own projects.

Notes:  This project was mainly developed to integrate ODEs along paths in complex domain, however, it can be simply modified to integrate ODEs along paths in (R^n), n >0, i.e., to integrate real dynamical systems.

Library

This section contains functions that performs numerical integration and plotting of solutions of ODEs (Ordinary Differential Equation s) along paths in some domain in comlex plane.
-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

Creates the list of equations along path form list of equations in above form

MakeSystem[list_,path_]:=Module[{equations,system},
equations=MapThread[D[#1[t],t]==D[path[t],t]*#2&,{list[[1]],list[[2]]}];
system=Flatten[{equations,list[[3]],{x'[t]==D[path[t],t],x[0]==path[0]}}];
system=Append[{system},Append[list[[1]],x]];
Return[system];
];

AreaQ[x_]:=Abs[x]≤5;  (*Area predicate*)
MaxNorm=1000; (*Maximal norm to stop integration*)

Path types (self explaining):

linePath[t_,pathbegin_,angle_]:=Exp[I*angle]*(t+pathbegin);
spiralPath[t_,pathbegin_, a_,b_,phi_,dir_]:=(pathbegin+(a*t+b)*Exp[I*dir*t])*Exp[I*phi]

Integrate along path with varable t along path, equation, area predicate and max norm of solution (Max(Abs[y],Abs[y’]] )  :

SolveAlongPath[t,equations_,path_,AreaQ_,MaxNorm_]:=Module[{pathSolution,sol,system},
system=MakeSystem[equations,path];
(*Print[system];*)
sol=NDSolve[system[[1]],system[[2]],{t,0,∞},Method->{"EventLocator","Event"->( Max[Abs[y'[t]],Abs[y[t]]]-MaxNorm>0|| ! AreaQ[x[t]]),"EventAction":>Throw[end=t,"StopIntegration"],"EventLocationMethod"->"LinearInterpolation",Method->"ExplicitRungeKutta"}];
Return[{{Re[x[t]],Im[x[t]],Abs[y[t]],Arg[y[t]]},end}/.sol[[1]]];
];

Plot solution along path as a paramteic curve colored by the phase function Hue[  (phase+π)/( 2π) ], where phase is in [-π,π]

PlotPathSolution[PathSolution_,AbsAxisName_:"AbsSol"]:=ParametricPlot3D[PathSolution[t][[1]][[1;;3]],{t,0,PathSolution[t][[2]]},AxesLabel->{Re,Im,AbsAxisName},ColorFunction->Function[t,Hue[(PathSolution[t][[1]][[4]]+π)/(2π)]] ]

Creates domain solution along paths emanating form the beginning point of the path as a function path[t,phi], where t is an parameter and parametrized by the angles form the range [phiStart, phisStop] with step of phiStep magnitude

DomainSolution[phiStart_,phiEnd_,phiStep_, t,equations_,path_,AreaQ_,MaxNorm_]:=Module[{tmpPath,domainSolution,PathSol,phi},
domainSolution={};
For[phi=phiStart,phi≤phiEnd,phi=phi+phiStep,tmpPath[t_]=path[t,phi];PathSol=SolveAlongPath[t,eq,tmpPath,AreaQ,MaxNorm];domainSolution=Append[domainSolution,PathSol];];
Return[domainSolution];
];

Show domain solution (sum of all paths) in new window to prevent gathering too many data in one file (problems with reading/writing large files).

(*ShowDomainSolution[domainSolution_]:=Module[{nb},
nb=CreateWindow[DocumentNotebook[{Show[Table[PlotPathSolution[#]&@Function[t,domainSolution[[i]]],{i,Length[domainSolution]}],PlotRange->All]}]];
Return[nb];
];*)

ShowDomainSolution[domainSolution_,Range_:All,AbsAxisName_:"AbsSol"]:=Module[{nb},
nb=CreateWindow[DocumentNotebook[{Show[Table[PlotPathSolution[#]&@Function[t,domainSolution[[i]]],{i,Length[domainSolution]}],PlotRange->Range]}]];
Return[nb];
];

Examples

Before you run this section execute Library section. It is recoomended to restart Kernel (Evaluation->QuitKernel or Quit[]), however you should first save you previous temporaty resluts as they will be lost when Kernel will be off.

Harmonic oscilator:

Define equation and initial data y’’[x]+y[x]==0 y[0]=1, y’[0]==0.   We integrate along path x[t], so the functions will depend on t not on x. Assume that you will pick the path that starts from x[0]==0. Then second order ODE can be rewritten as a system
y’[t]==v[t]
v’[t]==-y[t]
y[0]==1
v[0]==0
Write it in the form { { orderet dependent variables},{ordered in the same way rhs of equations}, {initial conditions} }

Warning:  It is assumed that y is a solution whoch you are interested to plot.

eq={{y,v},{v[t],-y[t]},{y[0]==1,v[0]==0}}

You can chceck how your system will be looking when written along path of integration

MakeSystem[eq,path]

Define the area of integration and maximum norm of solution (Max(Abs[y],Abs[y’]] ) which should not be crossed.

AreaQ[x_]:=Abs[x]≤5;  (*Area predicate*)
MaxNorm=1000; (*Maximal norm to stop integration*)

Define path:

path[t_]:=linePath[t,0,0];

Check that the path starts at initial data:

path[0]==0

Solve equation along path:

PathSolution[t_]=SolveAlongPath[t,eq,path,AreaQ,MaxNorm];

Plot a solution along path:

PlotPathSolution[PathSolution,"Abs"]

Path solution is a list  { { Re[ x[t] ],  Im[ x[t] ], Abs[ y[t] ], Arg[ y[t] ] }, maximal value of t}.   Maximal value of t indicates where integrations stops. It can be the point of boundary or the point where the solution crossed the norm constraint.

PathSolution[t]

You can now construct solution along many paths by filling a list of solutions using line path:

domainSolution={};
For[phi=0,phi<=2π,phi=phi+0.1,path[t_]=linePath[t,0,phi];PathSol=SolveAlongPath[t,eq,path,AreaQ,MaxNorm];domainSolution=Append[domainSolution,PathSol];]

or spiral path:

domainSolution={};
For[phi=0,phi<=2π,phi=phi+0.5,path[t_]=spiralPath[t,0,1,0,phi,1];PathSol=SolveAlongPath[t,eq,path,AreaQ,MaxNorm];domainSolution=Append[domainSolution,PathSol];]

and visualize it (it is commented because it is time consuming and generates large notebook file, which can crash Mathematica during saving this file):

(*Show[Table[PlotPathSolution[#]&@Function[t,domainSolution[[i]]],{i,Length[domainSolution]}],PlotRange->All]*)

It is mor efficient to create large plots in different notebook:

nb=CreateWindow[DocumentNotebook[{Show[Table[PlotPathSolution[#]&@Function[t,domainSolution[[i]]],{i,Length[domainSolution]}],PlotRange->All]}]];

It can be automatically done by the following three commands:

Create a path parametrized by “angle” (iteration variable along paths):

path[t_,phi_]:=linePath[t,0,phi];

Create domain solution in the range phi in [0,2π] with step 0.05

domainSolution=DomainSolution[0,2π,0.05, t,eq,path,AreaQ,MaxNorm];

Plot solution in new notebook:

ShowDomainSolution[domainSolution];

More elaborated example - the Lane-Emden equation

The equation is of the form y''[x]+(α/x)y[x]+((y[x])^p)==0, y[0]==1, y’[0]==0. We set α=2, p=5.   
x=0 is a fixed singuar point of equation coefficeint ((α/x))therefore to start integrate we have to shift initial data form fixed singularity using analytic solution derived below:

truncatedSeries=1+Sum[(Subscript[a, i])*(x^i),{i,2,6}];
subst=y''[x]+(2/x)*y'[x]+y'[x]+((y[x])^5)/.{y[x]->truncatedSeries,y'[x]->D[truncatedSeries,x], y''[x]->D[truncatedSeries,x,x]};
coeffList=CoefficientList[subst,x]//Simplify;
a2=Solve[coeffList[[1]]==0,(Subscript[a, 2])]
a3=Solve[(coeffList[[2]]/.a2)==0,(Subscript[a, 3])]

....  we are interested in illustrations therefore we limit ourselves only to the truncated expansion of the form 
y[x]=1-(1/6)(x^2), y'[x]=((-1)/3)x
The proof of the convergence of the series can be found in extensive literature of the subject.

We can now rewrite the equation as a system of first order ODEs.  Itroduce a path x[t], and treat y[x[t]] as simply y[t]:
y’[t]==v[t]
v’[t]==((-2)/(x[t]))v[t]-((y[t])^5)
y[0]==1-(1/6)((x[0])^2)
v[0]==((-1)/3)*x[0]

where x[0] is a beginning of the path in the radius of convergence of the power series (close to x==0).

eq={{y,v},{v[t],((-2)/(x[t]))v[t]-((y[t])^5)},{y[0]==1-(1/6)((x[0])^2),v[0]==((-1)/3)*x[0]}}

We take a path close to x==0.   In this setup x[0]==0.001*Exp[i * phase], e.g., |x[0]| == 0.001. Error of approximation of initial data is of order |x[0]((|)^4).

path[t_,phi_]:=linePath[t,0.001,phi]

Define domain border and maximum of norm

AreaQ[x_]:=Abs[x]≤5;  
MaxNorm=100;

Integrate along paths

domainSolution=DomainSolution[0,2π,0.1, t,eq,path,AreaQ,MaxNorm];

By simply plotting the domain solution we would get unreadable plot, becasue there are directions which “hit” or slightly miss singularites, e.g., path no.17:

GraphicsGrid[{Table[PlotPathSolution[#]&@Function[t,domainSolution[[i]]],{i,16,18}]}]

Therefore we should restrict the range of AbsSol to some small interval. This produce more readable output:

ShowDomainSolution[domainSolution,{All,All,{0,10}}];

Your project