Natural Solutions to Delay Differential Equations

James F Sturnfield, Consulting Mathematician

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

ABSTRACT

Mathematica has provided a strong platform for developing a new approach to solving problems containing delay differential equations (DDE). In particular, a method has been developed for finding “natural” analytic solutions that do not require a priori history to be specified. 

The “natural” solution will, in some sense, be determined from minimum energy requirements and from the DDE rather than an arbitrary choice of history. There will be some DDE where there can be multiple natural solutions, and additional conditions may be needed to determine an appropriate solution(s). It is also possible to determine natural solutions that satisfy a finite number of starting conditions. All other solutions to the DDE can be approximated by natural solutions corresponding to appropriate choices of starting conditions.

This approach does not require an arbitrary choice for a history and does not require as much data to distinguish between different models. An additional advantage of natural analytic solution usage is that solutions can be back calculated to an earlier time than the primary data set. This could help determine the start of specific features in a data set that might diagnose process problems or further validate the model. 

This presentation will illustrate the calculation of natural solutions for DDE with a single time delay for both linear and nonlinear equations with and without starting conditions being specified. In addition, examples of using this approach for fitting data of processes that include time delays will be demonstrated. In particular, some data from the epidemiology of COVID-19 will be examined.

Background

About Myself (James F Sturnfield)

James worked for Dow Chemical for 36 years in a wide range of various types of modeling. He is retired and working on projects of his own interest.

more information ...

He has a PhD from Purdue University in mathematics with a thesis that studied high energy singularities in a version of quantum field theory. In addition to studying functional differential equations, he focused on studying human movement and the impact of fascia layers. He has been using Mathematica for 21 years.

How I Came to this project

When Covid-19 was threatening lockdown in the United States, some of my conversations turned towards better understanding its spread.

Exchanges with Mark Kotanchek (Evolved Analytics https://evolved-analytics.com )

Notebooks by Robert B. Nachbar “Epidemiological Models for Influenza and COVID-19” [RB1] [RB4]

Delayed Differential Equations (DDE) seemed like a natural fit for modelling the transmission of COVID-19 because of the possibly long delays between the initial contagiousness and the symptoms/detection of the disease.
There seemed to be numerous challenges using these DDE in this problem.

Overview of Presentation

Background for delayed differential equations (DDE)

Explicit analytic solutions for some DDE

Introducing concept of natural solution to DDE

Original research development
	New terms used (concepts may have other other definitions or names)
		Natural Solution of DDE - “Analytic”  solution with minimal energy use 
		DDE minimal energy- Defined with appropriate metric (or rapid local convergence of Taylor series)
		Taylor Spline Fit - Polynomial spline curve that transitions between points with finite Taylor series
		DDE Keypoints - Points linked to initial point by DDE and time delay structure
		Linked Points (for single delay) - A given point with its lag point and lead point

Advantages of natural solutions

No history needed, just initial conditions
	No Discontinuities
	Back Calculation can be done
	Minimizes variations in the solution
	Possibly more realistic because of minimal energy input
	Improved Data Fitting of DDE

Tools towards finding natural solutions

Show how natural solutions might be used

Summary

Generate interest in further developing Natural Solutions

Quick Background of DDE

Delay differential equations (DDE) are differential types of equations where some of the functional terms in the equations involve argument at earlier times (smaller independent values). These equations can be written as a calculation of the derivative of a functional (unknown function) as a function of the time (independent variable), the functional at a time and the functional at various earlier times.

These equations are important in many different areas of engineering and science including viscoelasticity, nuclear reactors, microbiology, epidemiology, and gene transcription. These equations probably go back almost to the creation of calculus. The extensive studies of them seem to have begun in the late 1940's with the introduction of computers [HJ][RR][BC].

We will consider delay differential equations (DDE) that can be written in the form:
 g'[t]= 𝒢[t,g[t],g[t-τ]] where there is a single time delay τ.
When it is understood that g is a function, it will also be written as 
g'= 𝒢[t,g,(Subscript[g, τ])] where  (Subscript[g, τ])[t] is the function g[t -τ]. 
These can be easily extended to vector equations

Manipulate[
Module[{T = 3,td=1,a=8/9,sol,x,t,esp=0.00001},sol = First[x/. NDSolve[{x'[t] == a x[t](1 - x[t - td]), x[t /; t ≤ 0] == (3/4)+b t},x,{t,-td-1,T+1}]];
GraphicsRow[{Plot[sol[t],{t,-td-1,T}, PlotRange->{{-td-0.5,T},{-1,2}},PlotLabel-> "x[t]"],Plot[sol'[t],{t,-td-1,T}, PlotRange->{{-td-0.5,T},{-2,2}},PlotLabel-> "x'[t]",Exclusions->"Discontinuities",ExclusionsStyle->Blue],Plot[((sol'[t+esp]-sol'[t])/esp),{t,-td-0.5,T}, PlotRange->{{-td-1,T},{-3,2}},PlotLabel-> "x''[t]"]}]],{{b,0.5,"History Slope, b"},-1.5, 1.5,.25, Appearance-> "Open"},ControlPlacement->Bottom ,FrameLabel->Style["x'[t]= (8/9) x[t] (1 - x[t-1]) & for t<0, x[t]=(3/4)+b t",16,Bold] ]
closeCellGroup[]

DDE used in presentation

The linear DDE  (analytic solutions)

These are delay differential equations (DDE) of the form:
 g'[t]= α g[t]+ β g[t-τ]] where there is a single time delay τ.

Population growth is a non-linear DDE

The rate of growth can be modeled by k[𝒫]= (Subscript[k, 1])( ((M - 𝒫)/M)) where 𝒫 is the population size (a function of time), M is called the carrying capacity and (Subscript[p, 0 ])is the initial population. The time delay of the initial time is assumed to be within the known history of the population. The equation can be written as 𝒫'[t]=(Subscript[k, 1])(((M-𝒫[t]))/M)𝒫[t-τ].

Traditional approach to solving DDE

Utilizing ODE techniques

DDE are similar to an ordinary differential equation (ODE) and can almost be solved with similar approaches. The challenge is for a given point at an initial time (Subscript[t, 0]), the derivative at (Subscript[t, 0]) can not be calculated without the value at the delay time (Subscript[t, 0])- τ and these delay values are needed for each step of the integration process.

Traditional approach requires “history”

The traditional approach obtains a solution by requiring a “history” function which is given for the function behavior before the initial time. The solution is obtained with forward integration over progressive intervals of length τ that are solved sequentially.  The history replaces the delay functional for the integration interval. The solution for one interval becomes the history for the next interval.

For g'[t]= 𝒢[t,g[t],g[t-τ]] with g[t-τ]=ϕ[t] for (Subscript[t, 0])≤t≤(Subscript[t, 0])+τ,
the solution can be calculated for the first interval from the ODE g'[t]= 𝒢[t,g[t],ϕ[t]]with g[0]=ϕ[τ]

Example of traditional approach for population growth

Consider  𝒫'[t]=0.6(((250-𝒫[t]))/250)𝒫[t-2]with 𝒫[t-2]=ϕ[t] (ϕ[t]=100 for 0≤t≤2) as the “history” function, 
then 𝒫'[t]=0.6(((250-𝒫[t]))/250)ϕ[t]=0.6(((250-𝒫[t]))/250) 100 for  0≤ t ≤2 which can be solved as an ODE:

DSolveValue[{𝒫'[t]==60(1 - 𝒫[t]/250),𝒫[0]==100},𝒫,t][[2]]
closeCellGroup[]

This solution can be used for ϕ[t] for 2≤t≤4 to give an ODE that can find a solution for this range and this process is continued.

Comparison of different histories on traditional solutions

The choice of a priori history for finding a solution for a DDE will usually lead to discontinuities in the first and higher derivatives.

pExplicitFun=First[DSolve[{𝒫'[t]==(Subscript[k, 1]) 𝒫[t](1 - 𝒫[t]/M),𝒫[0]==(Subscript[p, 0])},𝒫,t, Assumptions-> 0< (Subscript[p, 0])&&0< M&&0<(Subscript[k, 1])]]; pExample[t_]:=(𝒫[t]/. pExplicitFun /. {(Subscript[p, 0])-> 100, M-> 250, (Subscript[k, 1])-> 0.6});
pDelayFix=First[NDSolve[{𝒫'[t]==(Subscript[k, 1]) 𝒫[t-τ](1 - 𝒫[t]/M),𝒫[t/;t<0]==(Subscript[p, 0])}/. {(Subscript[p, 0])-> 100, M-> 250, (Subscript[k, 1])-> 0.6,τ-> 2},𝒫,{t,-4,10}]];
pDelayAnalyticHistory=First[NDSolve[{𝒫'[t]==(Subscript[k, 1]) 𝒫[t-τ](1 - 𝒫[t]/M),𝒫[t/;t<0]==pExample[t]}/. {(Subscript[p, 0])-> 100, M-> 250, (Subscript[k, 1])-> 0.6,τ-> 2},𝒫,{t,-4,10}]];
GraphicsRow[{Plot[{pExample[t],𝒫[t]/.pDelayFix,𝒫[t]/.pDelayAnalyticHistory}, {t,-4,10},PlotLabel->"Population per Year",PlotStyle-> {{Red,Dashed},Blue,Green},PlotLegends->Placed[{"No time delay:ODE","2 Year delay:Fixed history","2 Year delay:ODE history"},Below],ImageSize->Medium],Plot[{𝒫'[t]/.pDelayFix,(((𝒫'[t+ϵ]-𝒫'[t]))/ϵ)/.pDelayFix/.{ϵ-> 0.000000001},
(* Note solution from NDSolve for a DDE has issues with the second derivative of the fixed history below initial time *)
Piecewise[{{𝒫'[t]/.pDelayAnalyticHistory,t≥ 0},{pExample'[t],t<0}}] ,(* Note NDSolve for DDE the derivative of the solution for this case does not equal the derivative of the history below the initial time *)Piecewise[{{(((𝒫'[t+ϵ]-𝒫'[t]))/ϵ)/.pDelayAnalyticHistory/.{ϵ-> 0.000000001},t≥ 0},{pExample''[t],t<0}}]},{t,-4,10},PlotLabel->"Derivatives of Population 𝒫[t] from different histories",PlotLegends->Placed[{"(\[DifferentialD]/(\[DifferentialD]t))𝒫[t] fixed history","((\[DifferentialD]^2)/(\[DifferentialD](t^2)))𝒫[t] fixed history","(\[DifferentialD]/(\[DifferentialD]t))𝒫[t] ODE history","((\[DifferentialD]^2)/(\[DifferentialD](t^2)))𝒫[t] ODE history"},Below],PlotStyle-> {Blue,{Blue,Dashed},Green,{Green,Dashed}},ImageSize->Medium]},Spacings->0]
closeCellGroup[2]

In the example above, the blue lines show the behavior of the solution with a fixed history function and the green lines show the behavior of the solution with the corresponding ODE being used for the history.

The integration process for solving the DDE smooths out the impact of history with each delay interval, so the impact of history is reduced with each interval step.

Sensitivity of different histories on solutions

For many DDE, changes in the history can create very different results to the solutions.

testF=Function[{a},NDSolveValue[Evaluate[{f'[t]== α f[t - δ]+β f[t],f[t/;t<=0]==1+a t}/.{α-> -15,δ->3,β->-5}],f,{t,-3,20}]];Plot[Evaluate[{testF[0][t],testF[1][t],testF[-1][t]}],{t,-3,6.5},PlotStyle->{Blue,{Red,Dashed},{Purple,Dashed}},PlotLegends->{"f[t<0]=1","f[t<0]=1+t","f[t<0]=1-t"},PlotLabel->"History impact on solutions for f'[t] = -5 f[t] -15 f[t-3]",ImageSize-> Large]
closeCellGroup[]

In the above plot, three different lines are used for the history with very different results.

Since arbitrary history can be specified for a solution to a DDE with a specified finite set of starting conditions, there are a continuum of solutions to  DDE.  Even if restricted to analytic solutions, there are a continuum of analytic solutions.

Note: There are also many DDE problems with given initial conditions, delay, and range of integration where the history has very little impact (this may be why this area of study has not been pursued).

Explicit analytic solutions for some DDE

Thanks to work of Clement Fabio [FC] and various students of A.G.Ulsoy [AU] (F.M.Asi)

Consider DDE of the form of f'[t]= α f[t-δ]+β f[t]  and solutions of the form 𝒻[t]=c (e^(d t)).  
These solutions can be shown to be 𝒻[t]= c (e^(((β+((ProductLog[(C[1]),(e^(-β δ)) α δ])/δ))t))).

Example with β=0

The ProductLog function is also known as the Lambert W function. It has multiple branches that give solutions in complex space. The (Subscript[c, 1]) indicates which branch is being considered. The principle branch is given by (Subscript[c, 1])=0. Because the DDE is linear and homogeneous, linear combinations of these solutions are also analytic solutions and both the real and imaginary parts of these solutions are also solutions.

A look at the general solutions for f'[t]= α f[t-δ]+β f[t]

Manipulate[
Module[{td=-1,T = 3,c=1,t},analyticFunSpec=c Re[(E^(t (β+((ProductLog[0,(E^(-β δ))α δ])/δ))))]; analyticFunSpecIm=c Im[(E^(t (β+((ProductLog[0,(E^(-β δ))α δ])/δ))))];
 analyticFun2Spec=c Re[(E^(t (β+((ProductLog[-1,(E^(-β δ))α δ])/δ))))];
 analyticFun2SpecIm=c Im[(E^(t (β+((ProductLog[-1,(E^(-β δ))α δ])/δ))))];Plot[{analyticFunSpec,analyticFunSpecIm,analyticFun2Spec,analyticFun2SpecIm},{t,td,T},PlotRange->{{td,T}(*,{-15,15}*)},PlotLabel-> "𝒻[t]= Exp \!\(TraditionalForm\`\([t(β + \*FractionBox[\(\(ProductLog[\\ α\\ δ\\ ]\)\(\\ \)\), \(δ\)])]\)\)",PlotStyle-> {Red,{Purple,Dashed},Blue,{Green,Dashed}},PlotLegends->{"𝒻[t] real","𝒻[t] Im","𝒻[t] 2nd branch real","𝒻[t] 2nd branch Im"}]]
,{{α,1,"Transition Rate, α"},-.2,3,.08, Appearance-> "Open"},{{β,1,"Growth Rate, β"},-3,2,.2, Appearance-> "Open"},{{δ,1.5,"Time delay, δ"},0.5,2,.01, Appearance-> "Open"},ControlPlacement->Bottom ,FrameLabel->Style["𝒻'[t]= α 𝒻[t-δ]+β f[t]",16,Bold] ]
closeCellGroup[]

The plot above looks at the real and imaginary parts of the principle branch and the first lower branch of the analytics solutions to the linear DDE. There are regions where there are two real solutions, only one real solution, and there are two complex solutions that are conjugates of each other.

Plot[Evaluate[{Exp[(β+((Re[ProductLog[(E^(-β δ)) α δ]])/δ))t]Cos[((Im[ProductLog[(E^(-β δ)) α δ]])/δ)t],Exp[(β+((Re[ProductLog[1,(E^(-β δ)) α δ]])/δ))t]Cos[((Im[ProductLog[1,(E^(-β δ)) α δ]])/δ)t],Exp[(β+((Re[ProductLog[2,(E^(-β δ)) α δ]])/δ))t]Cos[((Im[ProductLog[2,(E^(-β δ)) α δ]])/δ)t]}/.{α-> -8,δ->2,β->-6}],{t,-2.1,2.1},PlotLabel->"Analytic Solution to f'[t]= -8 f[t - 2]-6 f[t]" ]
closeCellGroup[]

The plot above shows three different analytic solutions to the DDE f'[t]= -8 f[t - 2]-6 f[t] with initial condition f[0]=1. The principle branch (blue curve) has the least variation of these cases and it is desirable to call it a more “natural” solution.

Introducing concept of natural solution to DDE

My original research development
	New terms introduced (concepts may have other definitions or names)

Natural solution of DDE - "analytic" solution with minimal “energy” use 
	DDE minimal energy- Defined with appropriate metric (or rapid local convergence of Taylor series)
	For ODE, a minimum curve can be defined by the metric determined by the ODE
	For DDE, the usual metric space is infinite dimension 
		○ Essentially because of the continuum of histories that can be defined
		○ Creates challenges for comparing two different solution curves
	Using concept of rate of convergence of Taylor series to nearby points
		○ Instead of using a metric, minimize local contribution of higher order terms in Taylor series.
		○ Will vary in different parts of range
		○ Still developing mathematical rigor
		○ Postulate that there is a corresponding metric

Also creating and utilizing new terms and concepts
	Taylor Spline - Polynomial spline curve that transitions between points with finite Taylor series
	Keypoints - Points linked to initial point by DDE and time delay structure
	Linked Points  - A given point with its lag point and lead point (for single delay)

Advantages of Natural Solutions to DDE

★ No history needed, just initial conditions
	○ There are also fundamental solutions that don’t need initial conditions
● No discontinuities 
	○ Approximations can be generated to any degree of continuity
● Back calculation can be done
	○ For many DDE, behavior of the solution can be determined for earlier times
● Minimizes variations in the solution
	○ Requirement of minimal energy reduces the possible variations
	○ Oscillations in solutions will occur in the lowest possible state
● Possibly more realistic because of minimal energy input
	○ Many physical systems tend to seek states of lowest energy requirements
● Improved data fitting of DDE
	○ As example, fitting “ history” data becomes an integrated part of obtaining the natural solution

Natural Analytic Solution for DDE

In determining what is a natural analytic solution for a DDE, there needs to be a “natural” way to build Taylor series at the initial, lag, and lead points (keypoints). In addition, there needs to be a way to build the curves between these points so that the integration of the DDE along these curves match the values at the lead Taylor series.

For example, if values are known for a set of keypoints, then DDE will give the derivative at the initial and the lead points and similarly the second derivative can be calculated at the lead point. The DDE solution can be approximated using curves between the key points that satisfy their derivatives.

Example: f'[t]= -8 f[t - 2]-6 f[t]  with f[0] =1. (f[-2]≈-0.7650242973906217`, f[2]≈-1.2372978296502908`)

(* General Linear DDE *)
ddeLEquation=D[f[t],t]==α f[t-τ] + β f[t];
tinitial=0;
(* Example equation f'[t]= -8 f[t - 2]-6 f[t] *)
eqLConstants={α-> -8,β-> -6, τ-> 2};
ddeLEquations=ddeLEquation/.eqLConstants;tlag=tinitial-τ/.eqLConstants;tlead=tinitial+τ/.eqLConstants;
(* Principle Branch *)
exampleSol=(E^(0.12019437912496667` t)) Cos[1.4541590317395714` t];

(* Analytic Solution to example *)
exampleSol=(E^(0.12019437912496667` t)) Cos[1.4541590317395714` t];pointList={-2,0,2};npoint=3;nDeriv={0,1,2};
(* get keypoints with increasing derivatives *)simpleProbKey=Table[Prepend[Table[(D[exampleSol,{t,i}]/.t->pointList[[j]]),{i,0,nDeriv[[j]]}],pointList[[j]]],{j,1,npoint}];
(* Calc Local Taylor series for keypoints to given derivatives *)
simpleTaylorApprox=Piecewise[Append[Map[{Sum[Rest[#][[i]] ((((t-First[#]))^(i-1))/((i-1)!)),{i,1,Length[Rest[#]]}],First[#]-ϵ<t≤First[#]+ϵ}&,simpleProbKey],{Missing,True}]];
(* generate polynomial spline curves continuous to degree at initial point *)
myCurve=Piecewise[{{a+ b t +c/2 (t^2)/.First[Solve[{a == simpleProbKey[[2,2]],b==simpleProbKey[[2,3]],a + b simpleProbKey[[1,1]] + (c/2) ((simpleProbKey[[1,1]] )^2)==simpleProbKey[[1,2]] },{a,b,c}]],simpleProbKey[[1,1]]<= t<= simpleProbKey[[2,1]]},{a+ b t +(c/2) (t^2)+(d/6)  (t^3)+(e/24)  (t^4)/.First[Solve[{a == simpleProbKey[[2,2]],b==simpleProbKey[[2,3]],a + b simpleProbKey[[3,1]] + (c/2) ((simpleProbKey[[3,1]] )^2)+ (d/6) ((simpleProbKey[[3,1]] )^3)+ (e/24) ((simpleProbKey[[3,1]] )^4)==simpleProbKey[[3,2]], b + c simpleProbKey[[3,1]] + (d/2) ((simpleProbKey[[3,1]] )^2)+ (e/6) ((simpleProbKey[[3,1]] )^3)==simpleProbKey[[3,3]] ,  c + d simpleProbKey[[3,1]] + (e/2) ((simpleProbKey[[3,1]] )^2)==simpleProbKey[[3,4]]  },{a,b,c,d,e}]],simpleProbKey[[2,1]]<= t<= simpleProbKey[[3,1]]}}];
(* Display *)
Show[Plot[simpleTaylorApprox/.{ϵ-> 1/4},{t,-2.25,2.25},PlotStyle->{Blue,Dashed},PlotLabel->"KeyPoints for f'[t]= -8 f[t - 2]-6 f[t]" ,PlotLegends-> {"Local Taylor Series"}],ListPlot[Evaluate[MapThread[Labeled[Take[#1, 2], #2] &, {simpleProbKey, {"Lag", "Init", "Lead"}}]],PlotLegends-> {"Key Point"},PlotStyle->{Red}],Plot[myCurve,{t,-2,2},PlotStyle->{Green},PlotLegends-> {"Polynomial Spline Fit of Key Point"}]]
closeCellGroup[]

The plot shows the original keypoints and the local Taylor series generated by the DDE. These local Taylor series are then connected by appropriate polynomials on each side of the initial point. The curve to the left of the initial point is the lag curve, and the curve to the right is the lead curve. These can be considered an initial approximation of the natural solution to the example DDE with its initial condition.

The forward integration of the DDE with the keypoint lag curve used as the “history” generates a “lag” DDE solution. It does not match the keypoint lead curve because the corresponding linked points (separated by the time delay) on the keypoint lag and lead curves do not satisfy the DDE.  The “lag” DDE solution does not depend on the keypoint lead curve in any way because of the way that ODEs are solved. It is not surprising that the lag solution does not end at the position of the lead point (although it is close).

GraphicsRow[{Show[Plot[myCurve,{t,-2,2},PlotStyle->{Green},PlotLegends-> Placed[{"Polynomial Fit of Key Point"},Below]],Plot[myCheck,{t,0,2},PlotStyle->{Red},PlotLegends-> Placed[{"DDE continued from lag"},Below]],ImageSize-> Medium],Show[Plot[myCurve,{t,1.5,2},PlotStyle->{Green}],Plot[myCheck,{t,1.5,2},PlotStyle->{Red}],ImageSize-> Small]}]

A (higher) derivative term for the lag point can be calculated that will create an appropriate history polynomial for the correct integration value at the lead point. This is seen in the “purple” curve on the right below.

(* Extending derivatives at lag point *)newCurve=Piecewise[{{a+ b t +c/2 (t^2)+(d/6)(t^3)/.First[Solve[{a == simpleProbKey[[2,2]],b==simpleProbKey[[2,3]],a + b simpleProbKey[[1,1]] + (c/2) ((simpleProbKey[[1,1]] )^2)+ (d/6) ((simpleProbKey[[1,1]] )^3)==simpleProbKey[[1,2]] , b  + c simpleProbKey[[1,1]] + (d/2) ((simpleProbKey[[1,1]] )^2)==lagDeriv},{a,b,c,d}]],simpleProbKey[[1,1]]<= t<= simpleProbKey[[2,1]]},{a+ b t +(c/2) (t^2)+(d/6)  (t^3)+(e/24)  (t^4)/.First[Solve[{a == simpleProbKey[[2,2]],b==simpleProbKey[[2,3]],a + b simpleProbKey[[3,1]] + (c/2) ((simpleProbKey[[3,1]] )^2)+ (d/6) ((simpleProbKey[[3,1]] )^3)+ (e/24) ((simpleProbKey[[3,1]] )^4)==simpleProbKey[[3,2]], b + c simpleProbKey[[3,1]] + (d/2) ((simpleProbKey[[3,1]] )^2)+ (e/6) ((simpleProbKey[[3,1]] )^3)==simpleProbKey[[3,3]] ,  c + d simpleProbKey[[3,1]] + (e/2) ((simpleProbKey[[3,1]] )^2)==simpleProbKey[[3,4]]  },{a,b,c,d,e}]],simpleProbKey[[2,1]]<= t<= simpleProbKey[[3,1]]}}];
(* Calculate DDE solution for lead point *)
newCheck= f/.First[ParametricNDSolve[{ddeLEquations,f[t/;t<= 0]==newCurve},f,{t,0,2},lagDeriv]];
(* Adjust the lagderiv to reach lead point *)
newDeriv=FindRoot[newCheck[lagDeriv][2]-simpleProbKey[[3,2]],{lagDeriv,0}];
Show[Plot[{myCurve,newCurve/.newDeriv},{t,-2,2},PlotStyle->{Green,{Blue,Dashed}},PlotLegends->{Placed[ {"Polynomial Fit of Key Point"},Below], Placed[{"Adjusted Polynomial Fit"},Below]},PlotLabel-> "Lag pt deriv added"],
Plot[{myCheck,Evaluate[newCheck[lagDeriv][t]/.newDeriv]},{t,0,2},PlotStyle->{Red,Purple},PlotLegends->{ Placed[{"DDE check from orginal lag"},Below], Placed[{"DDE check from new lag"},Below]}],ImageSize->Small]

Show[Plot[myCurve-myCheck,{t,0,2},PlotStyle->{Red},PlotLegends-> Placed[{"Polynomial Fit - Original Check"},Above]],Plot[myCurve-Evaluate[newCheck[lagDeriv][t]/.newDeriv],{t,0,2},PlotStyle->{Purple},PlotLegends-> Placed[{"Polynomial Fit - newCheck"},Above],PlotLabel-> "lead fit - integration check"]]

The new DDE integration will then reach the lead point although there will still be some error in the shape of both the lag and lead curve.

Using “Taylor” splines to refine the curves

The polynomial curves used in the previous example will require higher orders of time to satisfy increasing derivatives at the key points. This increases the “energy” of the solution because higher order of Taylor series are being introduced for the points between these key points. In order to keep the approximation close to the “natural” solution, a spline approximation is introduced that will have polynomial pieces that are, at most, one order of derivative higher than the keypoints.

Polynomial spline knots can be inserted between each key position so that the spline satisfies the values and derivatives already determined at these key points and is continuous to at least the derivative order of the less defined key position.

plValues=Map[{Re[#],Im[#]}&,Map[N[ProductLog[(C[1]),(E^(-β δ)) α δ]]/.{α-> -8,δ->2,β->-6,(C[1])-> #}&,{0,1,2}]];exampleSolutions=Flatten[Map[{ Exp[(β+((#[[1]])/δ))t]Cos[((#[[2]])/δ)t],Exp[(β+((#[[1]])/δ))t]Sin[((#[[2]])/δ)t]}/.{α-> -8,δ->2,β->-6}&,plValues]];pointList={-2,0,2};npoint=3;nDeriv={1,2,3};linearProbKey=Map[Table[Prepend[Table[(D[#,{t,i}]/.t->pointList[[j]]),{i,0,nDeriv[[j]]}],pointList[[j]]],{j,1,npoint}]&,exampleSolutions];
linearTaylorApprox=Map[Piecewise[Append[#,{Missing,True}]]&,Map[{Sum[Rest[#][[i]] ((((t-First[#]))^(i-1))/((i-1)!)),{i,1,Length[Rest[#]]}],First[#]-ϵ<t≤First[#]+ϵ}&,linearProbKey,{2}]];
extendTSpt=Map[{#,linearTaylorApprox[[1]]/.{ϵ->1/2,t-> #}}&,{-5/3,-1/3,1/3,5/3}];splineApprox=Map[splineTaylorPoints[#,0]&,linearProbKey];splineProbKeyTemp=Map[{Prepend[Table[D[#[[1]],{t,i}](i!)/.{t->#[[2,1]]},{i,0,Length[#[[1]]]-1}],#[[2,1]]],Prepend[Table[D[#[[1]],{t,i}](i!)/.{t->#[[2,-1]]},{i,0,Length[#[[1]]]-1}],#[[2,-1]]]}&,Drop[First[Drop[Apply[List,splineApprox[[1,-1]]],-1]],1]];splineProbKey=Append[Map[Drop[First[#],-1]&,splineProbKeyTemp],Drop[Last[Last[splineProbKeyTemp]],-1]];splineTaylorApprox=Map[Piecewise[Append[#,{Missing,True}]]&,Map[{Sum[Rest[#][[i]] ((((t-First[#]))^(i-1))/((i-1)!)),{i,1,Length[Rest[#]]}],First[#]-ϵ<t≤First[#]+ϵ}&,{splineProbKey},{2}]];Show[{Plot[Evaluate[First[splineTaylorApprox]/.{ϵ-> 1/6}],{t,-2,2},PlotStyle->{Purple,Dashed},PlotLegends-> {"Local Taylor series"},PlotLabel->"Solving f'[t]= -8 f[t - 2]-6 f[t]"],ListPlot[Evaluate[MapThread[Labeled[Take[#1, 2], #2] &, {splineProbKey, {"Lag","Node","Node", "Init","Node","Node","Node", "Lead"}}]],PlotLegends-> {"Spline Points"}]}]

Show[Plot[splineApprox[[1,-1]],{t,-2,2},PlotStyle->Red,PlotLegends->{"Spline Approx"},PlotLabel->"Solving f'[t]= -8 f[t - 2]-6 f[t]" ],Plot[exampleSolutions[[1]],{t,-2,2},PlotStyle->Green,PlotLegends->{"Principal Solution"}],
Plot[Evaluate[First[splineTaylorApprox]/.{ϵ-> 1/4}],{t,-2-1/4,2+1/4},PlotStyle->{Purple,Dashed},PlotLegends->{"Spline node info"}]]

Note this process can be continued on each piece of the spline approximation and with higher order derivatives at each point (with higher order polynomials).

As the derivatives at the key point are adjusted, this “Taylor Spline” will be generated to produce the curve.

Back calculation

The forward integration of the lag approximation curve gives a DDE solution for that lag curve. In many cases, a back calculation can be done to extend the lead approximation curve into the lag interval for a lead DDE solution.
In general, if the value and its derivative are known at some position, then the DDE creates an implicit function for the value at the lag position point that could have one, multiple, or no solutions.

Example: Single solution 
		f'[t]= -8 f[t - 2]-6 f[t],   f[0]= 1 and f'[0]=0.120194,  then  f[t-2]= ((-6 f[t]-f'[t])/8)and f[-2]= -0.765024

Example: Multiple solutions  
		f'[t]= ((( f[t - 2]))^2)- f[t],   f[0]= 2 and f'[0]=-1,  then  f[t-2]= ±(Sqrt[f[t]+ f'[t]]) and f[-2]= ±1 ( note for f'[0]<-2, no real solutions)

Example: Singularity  
		𝒫'[t]=0.6(((250-𝒫[t]))/250)𝒫[t-2],   f[0]= 250 and f'[0]=-1,  then 𝒫[t-2]=(1250/3)((𝒫'[t])/((250-𝒫[t])))and 𝒫[-2][-2] is a singularity

The lead spline curve from the previous slide can be back calculated to give  a history that produces the lead curve with the DDE .

Extending the DDE condition to nearby point

Local keypoint behavior of DDE

DDE forces relationships between starting points and their lag points. There still needs to be additional relations from the DDE between the starting points and their lead points. This comes from the integration of the lag curve. Using the local Taylor Series at the key points, the values at nearby points can be estimated.

Show[Plot[exampleSolutions[[1]],{t,-2,2},PlotStyle->Yellow,PlotLegends->{"Principal Solution"},PlotLabel->"Solving f'[t]= -8 f[t - 2]-6 f[t]"],
Plot[Evaluate[linearTaylorApprox[[1]]/.{ϵ-> 1/3}],{t,-2,2},PlotStyle->{Purple,Dashed},PlotLegends->{"Current node info"},PlotLabel->"Solving f'[t]= -8 f[t - 2]-6 f[t]"],ListPlot[Evaluate[MapThread[Labeled[Take[#1, 2], #2] &, {simpleProbKey, {"Lag", "Init", "Lead"}}]],PlotLegends-> {"Key Point"},PlotStyle-> Red],ListPlot[Evaluate[MapThread[Labeled[#1, #2] &, {extendTSpt, {"Lag+", "Init-","Init+", "Lead-"}}]],PlotLegends-> {"Local Taylor step"},PlotStyle-> Blue]]
closeCellGroup[]

In particular, the DDE yields the derivative at the points near the initial position from the values for the solution near the lag point. The derivative at the points near the lead position can be calculated from the value near the initial point. If it is possible to back calculate the DDE, then similar calculations can determine similar information near the initial points from the values near the lead point. These calculation will not necessarily arrive at the same results near the initial points. Higher derivatives can be added at the lag point that will drive these competing methods to the same value.

linkedPoints={{-(5/3),-0.7075674004711744`,0.17237069075834188`},{-(1/3),0.8432612180240927`,0.8202383127304771`,-2.1001318008165315`},{(1/3),0.9233908041074038`,-0.5798495544805438`,-2.1001318008165315`},{(5/3),-0.9178515624859765`,-1.3066365103507178`,1.8868178299251133`,1.8267158089938493`}};
linkedArray=Map[MapAt[First[#]&,#,{1}]&,Map[{#}&,linkedPoints,{2}]];
{debug,ddeNearLag}=forwardCalcPtDDE[{ddeLEquations},Take[linkedArray,1],Map[Take[#,2]&,Take[linkedArray,{-2}]],{t,f}];
{debug,ddeNearInit}=forwardCalcPtDDE[{ddeLEquations},Take[linkedArray,{2}],Map[Take[#,2]&,Take[linkedArray,-1]],{t,f}];
ddeNear=Join[ddeNearLag,ddeNearInit];
taylorNearApprox=Map[Piecewise[Append[#,{Missing,True}]]&,Map[{Sum[Rest[#][[i]] ((((t-First[#]))^(i-1))/((i-1)!)),{i,1,Length[Rest[#]]}],First[#]-ϵ<t≤First[#]+ϵ}&,{ddeNear},{2}]];
Show[Plot[exampleSolutions[[1]],{t,-2,2},PlotStyle->Yellow,PlotLegends->{"Principal Solution"},PlotLabel->"Solving f'[t]= -8 f[t - 2]-6 f[t]"],
Plot[Evaluate[linearTaylorApprox[[1]]/.{ϵ-> 1/3}],{t,-2,2},PlotStyle->{Purple,Dashed},PlotLegends->{"Current node info"},PlotLabel->"Solving f'[t]= -8 f[t - 2]-6 f[t]"],ListPlot[Evaluate[MapThread[Labeled[Take[#1, 2], #2] &, {simpleProbKey, {"Lag", "Init", "Lead"}}]],PlotLegends-> {"Key Point"},PlotStyle-> Red],ListPlot[Evaluate[MapThread[Labeled[#1, #2] &, {extendTSpt, {"Lag+", "Init-","Init+", "Lead-"}}]],PlotLegends-> {"Local Taylor step"},PlotStyle-> Blue],Plot[Evaluate[taylorNearApprox/.{ϵ-> 1/4}],{t,-2,2},PlotStyle->{Green,Dashed},PlotLegends->{"DDE solve of near lead node"},PlotLabel->"Solving f'[t]= -8 f[t - 2]-6 f[t]"]]
closeCellGroup[]

Population DDE example for natural solution

A simple population DDE model will be used to illustrate the process of developing a natural solution.

where M=250 is the carrying capacity, (Subscript[k, 1])=0.6 is the maximum growth rate, and τ = 2 is the time delay and considered natural solutions without any additional conditions.

Initially, a quadratic approximation is made through the key points -2, 0, and 2. There are five possible natural solutions. Here the focus will only be on the two middle cases.

ddeEquation=D[𝒫[t],t]==(Subscript[k, 1])(((M-𝒫[t]))/M)𝒫[t-τ];

eqConstants={M-> 250,(Subscript[k, 1])-> 3/5, τ-> 2};
ddeEquations=ddeEquation/.eqConstants;tinitial=0;tlag=tinitial-τ/.eqConstants;tlead=tinitial+τ/.eqConstants;
approxForm= (Subscript[a, 0])+t (Subscript[a, 1])+(((t^2) (Subscript[a, 2]))/2);approxSub=Table[(D[𝒫[t],{t,i}]/.t-> t_)-> D[approxForm,{t,i}],{i,2,0,-1}];
approxEq=Flatten[Table[Table[D[ddeEquations,{t,j}]/.approxSub/.{t->i τ},{i,j,1}],{j,0,1}]/.eqConstants];
populationSolve=SortBy[Map[approxForm/.#&,NSolve[approxEq,{(Subscript[a, 0]),(Subscript[a, 1]),(Subscript[a, 2])}]],N[#]/.{t->0}&];
forwardIntChk=Map[Flatten[NDSolve[{ddeEquations,𝒫[t/;t<= tinitial]==#},𝒫, {t,tlag,tlead}]]&,populationSolve];GraphicsRow[{Show[Plot[Evaluate[Map[𝒫[t]/.#&, forwardIntChk[[2;;5]]]],{t,tlag,tlead},PlotStyle->{Red,Blue,Green,Purple},PlotLegends->Placed[{"𝒫[0]=0","𝒫[0]=50.026","𝒫[0]=196.137","𝒫[0]=250"},Below],PlotLabel->"Forward Integration: Pop DDE",PlotRange->{All,{-0.1,250.1}}],Plot[{populationSolve[[3]],populationSolve[[4]]},{t,0,tlead},PlotStyle->{{Cyan,Dashed},{Orange,Dashed}},PlotRange->{All,{-0.1,250.1}}]],
Plot[Evaluate[{((100(( 𝒫[t]/. forwardIntChk[[3]])- populationSolve[[3]]))/(populationSolve[[3]])),((100((𝒫[t]/. forwardIntChk[[4]])- populationSolve[[4]]))/(populationSolve[[4]]))}],{t,tinitial,tlead},PlotStyle->{Blue,Green},PlotLegends->Placed[{"𝒫[0]=50.026","𝒫[0]=196.137"},Below],PlotLabel->"% Error difference forward integration and quad approx: Pop DDE"]}]

The dashed lines in the first plot show the quadratic approximation while the solid lines are results from integrating these approximations for the lead part of a solution. The second plot shows the difference between the forward integration and the initial approximation.

Refining the approximation for the natural solution

The DDE can also be used as an implicit function in order to back calculate the lag section of a solution from the lead part of the approximate. Comparisons of forward integration and back calculation with the current approximation gives a measure of error on the fit with the DDE.

backSol=Map[Simplify[First[NSolve[{Evaluate[D[𝒫[t]/.{𝒫[t]->N[#]},t]==(ddeEquations[[2]]/.{𝒫[t+x_]-> (Subscript[𝒫, -x])})/.{𝒫[t]->N[#]}]},(Subscript[𝒫, 2])]]/.{(Subscript[𝒫, 2])-> 𝒫,t-> t+2}]&,populationSolve];
GraphicsRow[{Show[Plot[Evaluate[{𝒫/. backSol[[3]],𝒫/. backSol[[4]]}],{t,tlag,tinitial},PlotStyle->{Blue,Green},PlotLegends->Placed[{"𝒫[0]=50.026 back calc","𝒫[0]=196.137 back calc","𝒫[0]=50.026 quadratic","𝒫[0]=196.137 quadratic"},Below],PlotLabel->"Back Calculation: Pop DDE",ImageSize->Medium],Plot[Evaluate[{ populationSolve[[3]],populationSolve[[4]]}],{t,tlag,tlead},PlotStyle->{{Cyan,Dashed},{Orange,Dashed}},PlotLegends->Placed[{"𝒫[0]=50.026 quadratic","𝒫[0]=196.137 quadratic"},Below]], PlotRange->All,ImageSize->Medium]
,Plot[Evaluate[{100(((𝒫/.  backSol[[3]])- populationSolve[[3]])/(populationSolve[[3]])),100(((𝒫/.  backSol[[4]])-populationSolve[[4]])/(populationSolve[[4]]))}],{t,tlag-0.1,tinitial+0.1},PlotStyle->{Blue,Green},PlotLegends->Placed[{"𝒫[0]≈50 Calc - Quad","𝒫[0]≈196 Calc - Quad"},Below],PlotLabel->"% error in difference Back Calculation and Quad Approx: Pop DDE",ImageSize->Medium]}]

The solid line in the first plot shows the back calculated solution to the DDE. The second plot shows the errors between the initial approximation and the DDE calculations from that approximation. Note the difference at the lag and initial key points are zero while the difference at the lead point is caused by the shape of the connecting curve and the integration process.

In order to reduce this error, a higher derivative term is introduced at the lag point and extended points are introduced in the lag region, so the new approximation curve satisfies the DDE at the extended point and the integral between linked points in the lag section reach the corresponding lead point.

For the population models being considered, the adjustment to the original quadratic approximation is fairly small.

approxForm3= (Subscript[a, 0])+t (Subscript[a, 1])+(((t^2) (Subscript[a, 2]))/2)+(((t^3) (Subscript[a, 3]))/6);approxCond3=Transpose[Flatten[Table[Table[Map[D[approxForm3,{t,j}] ==#&, D[populationSolve,{t,j}]]/.{t->(i-1) τ},{i,j,1}],{j,0,1}]/.eqConstants,1]];approxEq3=Map[First[Solve[#,{(Subscript[a, 0]),(Subscript[a, 1]),(Subscript[a, 3])}]]&,approxCond3];forwardParmIChk=Map[Flatten[ParametricNDSolve[{ddeEquations,𝒫[t/;t<= tinitial]==(approxForm3/. #)},𝒫, {t,tinitial,tlead},{(Subscript[a, 2])}]]&,approxEq3];
approxEqA2=MapThread[(𝒫[(Subscript[a, 2])][2]/.#1)-(#2/.{t-> 2})&,{forwardParmIChk,populationSolve}]; 
(* Can't get solutions with NSolve or NSolveValues, so using FindRoot *)
Off[FindRoot::lstol];
{populationSolveLagNew,solveErr}= Transpose[MapThread[(rootFound=First[FindRoot[#2,{(Subscript[a, 2]),0}]];{Simplify[(approxForm3/.#1)/.rootFound],Simplify[#2/.rootFound]})&,{approxEq3,approxEqA2}] ];
On[FindRoot::lstol];
Show[Plot[Evaluate[populationSolveLagNew[[3;;4]]],{t,tlag,tinitial},PlotStyle->{Blue,Red},PlotLegends->{"𝒫[0]=50.026 quadratic","𝒫[0]=196.137 quadratic"},PlotLabel->"Refining Lag curve: Pop DDE"],Plot[Evaluate[populationSolve[[3;;4]]],{t,tlag,tinitial},PlotStyle->{{Cyan,Dashed},{Orange,Dashed}},PlotLegends->{"𝒫[0]=50.026 original","𝒫[0]=196.137 original"}], PlotRange->All]
closeCellGroup[]

The resulting error in the lead section is greatly reduced.

Show[Plot[Evaluate[{((100(( 𝒫[t]/. forwardParmIChk[[3]])- populationSolve[[3]]))/(populationSolve[[3]])),((100((𝒫[t]/. forwardParmIChk[[4]])- populationSolve[[4]]))/(populationSolve[[4]]))}],{t,tinitial,tlead},PlotStyle->{Blue,Green},PlotLegends->{"𝒫[0]=50.026 new lag curve","𝒫[0]=196.137 new lag curve"},PlotLabel->"% Error difference forward integration and quad approx: Pop DDE",PlotRange-> All],Plot[Evaluate[{((100(( 𝒫[t]/. forwardIntChk[[3]])- populationSolve[[3]]))/(populationSolve[[3]])),((100((𝒫[t]/. forwardIntChk[[4]])- populationSolve[[4]]))/(populationSolve[[4]]))}],{t,tinitial,tlead},PlotStyle->{{Blue,Dashed},{Green,Dashed}},PlotLegends->{"𝒫[0]=50.026 orginal lag","𝒫[0]=196.137 original lag"},PlotLabel->"% Error difference forward integration and quad approx: Pop DDE"]]
closeCellGroup[];

Natural solution for Delayed SIR Model

The SIR model of epidemiology put the population being considered into three different compartments: susceptible, infectious, and recovered. The rate of change in each compartment is related to the transference population between compartments. Time delays can occur in a variety of ways. In this example, we consider a delay effect on the transfer from the infectious group to the recovered group.

Equations for the delay SIR model for epidemiology

where the time delay τ = (1/4), β = 0.005028, and γ= 1.938. These values were calculated so that the initial behavior of this time delayed equation corresponded to the initial behavior of an ODE originally used to model the system.

The delay SIR model has significant differences when compared with a corresponding ODE model. In this case, there are not significant differences in the results for the natural solution and the use of a fixed history.

Fitting COVID-19 data with natural DDE models

Notes about the dataset

Using the Hubei, China data set provided through Wolfram, a fitting can be made using a model similar to the delay SIR model with some expanded compartments and parameters. There are a couple of corrections that were made to the data set which created jumps in the different quantities.

DateListPlot[{covidHubeiListex[[All,1;;2]],covidHubeiListex[[All,{1,3}]],covidHubeiListex[[All,{1,4}]],covidHubeiListex[[All,{1,5}]],covidHubeiListex[[All,{1,6}]]},PlotLabel->"Hubei",FrameLabel->{"time (d)","# individuals"},PlotLegends->{"Detected","Quarantined","Recovered","Dead","Removed"},Joined-> False,PlotRange->{{{2020,1,22},{2020,5,1}},Full},FrameTicks->{{Automatic,Automatic},{{{2020,2,1},{{2020,2,13},"Def change",.05,Red},{2020,3,1},{2020,4,1},{2020,5,1}}},None}]

Notes:
Removed = recovered + death
Quarantine = confirmed - recovered - death
(DateObject["Thu 23 Jan 2020", DateObject[{2020,1,23},"Day","Gregorian",-4.`]]) Wuhan in lockdown (beginning of this data set)
 (DateObject["Sat 15 Feb 2020", DateObject[{2020,2,15},"Day","Gregorian",-4.`]]) changed guidelines for determining COVID to include non DNA (large jump in confirmed cases)
  (DateObject["Sun 23 Feb 2020", DateObject[{2020,2,23},"Day","Gregorian",-4.`]]) A  variety of additional changes are made to reports over this weekend
(DateObject["Fri 17 Apr 2020", DateObject[{2020,4,17},"Day","Gregorian",-4.`]]) Additional deaths discovered and reported

Model and fitting

For this model, the compartments are: total detected cases, infectious cases (not measured because sign of COVID not seen at initial infection), quarantine cases(assuming cases are put into quarantine as soon as first detected), and removed cases (either recovered or dead). A time delay is included between first infectious until detected and quarantined.

The fit is a weighted sum of least squares based on expected quality of the data.

forceOfInfectionSIQR={λ[t_]:>((ℐ[t]+α 𝒬[t])/(𝒩-(1-α)𝒬[t]))};
wuhanSIQR={ℐ'[t]==-δ ℐ[t-τ]-μ ℐ[t]+β (𝒩 -ℐ[t]-𝒬[t]-ℛ[t]) λ[t],𝒬'[t]==δ ℐ[t-τ]-γ 𝒬[t],ℛ'[t]==γ 𝒬[t]+μ ℐ[t]}/.forceOfInfectionSIQR;
varsSIQR={ℐ,𝒬,ℛ};
infectInit=iDelay+(iInit-iDelay)(t+τ +24)/τ;wuhanParm={𝒩-> 80000,τ-> 6,α-> 0.05,δ-> .10,γ-> 0.05,μ-> 0.001,β-> .4};initCondWuhan0={ℐ[t/;t<= -24]==infectInit,𝒬[t/;t<= -24]==399,ℛ[t/;t<= -24]==45};
mySolveParm=ParametricNDSolve[Evaluate[Join[wuhanSIQR,initCondWuhan0]],varsSIQR,{t,-24,36},{𝒩,τ,α,δ,γ,μ,β,iInit,iDelay}];
{infectedSIQR,quarantinedSIQR,removedSIQR}={ℐ,𝒬,ℛ}/.mySolveParm;
detectedSIQR=Plus[infectedSIQR,quarantinedSIQR,removedSIQR];
susceptibleSIQR=𝒩-infectedSIQR-quarantinedSIQR-removedSIQR;
covidHubeiListCt=MapIndexed[Prepend[Drop[#1,1],First[#2]-25]&,covidHubeiListex];
Show[
Plot[Evaluate[{
detectedSIQR[pop,tau,alpha,delta,gamma,mu,beta,iInit,iDelay][t],
infectedSIQR[pop,tau,alpha,delta,gamma,mu,beta,iInit,iDelay][t],
quarantinedSIQR[pop,tau,alpha,delta,gamma,mu,beta,iInit,iDelay][t],
removedSIQR[pop,tau,alpha,delta,gamma,mu,beta,iInit,iDelay][t]
}/.{pop-> 80000,tau-> 6,alpha->0.0491298303917148`,delta->0.061409827798194226`,gamma->0.05440040263372297`,mu->0.001`,beta->0.4`,iInit-> 3400,iDelay-> 400}],{t,-24,36},
PlotRange->{{-24,36},{0,70000}},
PlotRangePadding->{Automatic,{Scaled[0.03],Scaled[0.03]}},
PlotLegends->{"Detect","ℐ","𝒬","ℛ"}(*,
PlotLabel->localDead[pop,10,beta,zeta,gamma,delta,mu][80]*),
FrameLabel->{"time","# people"},(*
GridLines->{None,(*{0}*)None},*)
PlotStyle->{
ColorData["Crayola"]["Green"],
ColorData["Crayola"]["Scarlet"],
ColorData["Crayola"]["Fuchsia"],
ColorData["Crayola"]["NavyBlue"]
}
],
ListPlot[Transpose[covidHubeiListCt[[1;;61,{2,3,6}]]],DataRange->{-24,36},
PlotStyle->{ColorData["Crayola"]["Green"],ColorData["Crayola"]["Fuchsia"],
ColorData["Crayola"]["NavyBlue"]}
]
]

The above model is not detailed enough to make a good fit of the data.

Show[
Plot[Evaluate[{
detectedSIQR[pop,tau,alpha,delta,gamma,mu,beta,iInit,iDelay][t],
infectedSIQR[pop,tau,alpha,delta,gamma,mu,beta,iInit,iDelay][t],
quarantinedSIQR[pop,tau,alpha,delta,gamma,mu,beta,iInit,iDelay][t],
removedSIQR[pop,tau,alpha,delta,gamma,mu,beta,iInit,iDelay][t]
}/.{pop-> 80000,tau-> 6,alpha-> 0.05,delta-> .10,gamma-> 0.05,mu-> 0.001,beta-> .4,iInit->3400,iDelay-> 400}],{t,-24,36},
PlotRange->{{-24,36},{0,70000}},
PlotRangePadding->{Automatic,{Scaled[0.03],Scaled[0.03]}},
PlotLegends->{"Detect","ℐ","𝒬","ℛ"}(*,
PlotLabel->localDead[pop,10,beta,zeta,gamma,delta,mu][80]*),
FrameLabel->{"time","# people"},(*
GridLines->{None,(*{0}*)None},*)
PlotStyle->{
ColorData["Crayola"]["Green"],
ColorData["Crayola"]["Scarlet"],
ColorData["Crayola"]["Fuchsia"],
ColorData["Crayola"]["NavyBlue"]
}
],
ListPlot[Transpose[covidHubeiListCt[[1;;61,{2,3,6}]]],DataRange->{-24,36},
PlotStyle->{ColorData["Crayola"]["Green"],ColorData["Crayola"]["Fuchsia"],
ColorData["Crayola"]["NavyBlue"]}
]
]

An alternative fit of the quarantine data does a better fit. In either model, the time delay equations do allow an estimation of the number of people who were infectious, but not yet detected.

DynamicModule[{localData,localSusceptible,localDetected,localInfected,localQuarantined,localRemoved,tau,alpha,delta,gamma,mu,beta,tMax=100,scale},
localData=Transpose[covidHubeiListCt[[1;;61,{2,3,6}]]];
{localSusceptible,localInfected,localQuarantined,localRemoved}={susceptibleSIQR,infectedSIQR,quarantinedSIQR,removedSIQR};
localDetected=Plus[localQuarantined,localRemoved];
Manipulate[
tau=inTau;alpha=10.^logAlpha;delta=10.^logDelta;
gamma=10.^logGamma;mu=10.^logMu;beta=10.^logBeta;
scale=50000/pop;
Show[
Plot[{
Plus[localQuarantined[pop,tau,alpha,delta,gamma,mu,beta,iInit,iDelay][t],
localRemoved[pop,tau,alpha,delta,gamma,mu,beta,iInit,iDelay][t]],
localInfected[pop,tau,alpha,delta,gamma,mu,beta,iInit,iDelay][t],
localQuarantined[pop,tau,alpha,delta,gamma,mu,beta,iInit,iDelay][t],
localRemoved[pop,tau,alpha,delta,gamma,mu,beta,iInit,iDelay][t]
},{t,-24,tmax},
PlotRange->{{-24,36},{0,ymax}},
PlotRangePadding->{Automatic,{Scaled[0.03],Scaled[0.03]}},
PlotLegends->{"Detect","ℐ","𝒬","ℛ"}(*,
PlotLabel->localDead[pop,10,beta,zeta,gamma,delta,mu][80]*),
FrameLabel->{"time","# people"},(*
GridLines->{None,(*{0}*)None},*)
PlotStyle->{
ColorData["Crayola"]["Green"],
ColorData["Crayola"]["Scarlet"],
ColorData["Crayola"]["Fuchsia"],
ColorData["Crayola"]["NavyBlue"]
}
],
ListPlot[localData,PlotMarkers->{Automatic,Small},
PlotStyle->{ColorData["Crayola"]["Green"],ColorData["Crayola"]["Fuchsia"],
ColorData["Crayola"]["NavyBlue"]},DataRange->{-24,36}
]
],
{{pop,80000,"effective population size"},1000{60,65,70,75,80,85,90,95,100,105}},
{{logBeta,Log10@.4,Dynamic@StringForm["β = `` (infection)",TraditionalForm[10.^logBeta]]},Log10@.01,Log10@5.0},
{{logMu,Log10@0.001,Dynamic@StringForm["μ = `` (infection remove)",TraditionalForm[10.^logMu]]},Log10@(0.0005),Log10@(0.1)},
{{logGamma,Log10@0.05,Dynamic@StringForm["γ = `` (remove rate)",TraditionalForm[10.^logGamma]]},Log10@(1/100.),Log10@(1/5.)},
{{logDelta,Log10@0.1,Dynamic@StringForm["δ = `` (quarantine rate)",TraditionalForm[10.^logDelta]]},Log10@(1/100.0),Log10@(1/0.05)},
{{logAlpha,Log10@0.05,Dynamic@StringForm["α = `` (quarantine infect)",TraditionalForm[10.^logAlpha]]},Log10@(1/1000),Log10@(1/10)},
{{inTau,6.,Dynamic@StringForm["τ = `` (detection delay)",TraditionalForm[inTau]]},3,20},
{{iInit,3400,Dynamic@StringForm["iInitial = `` (initial infection)",TraditionalForm[iInit]]},1000,6000},
{{iDelay,400.,Dynamic@StringForm["idelay = `` (delayed infection)",TraditionalForm[iDelay]]},0,1000},
Delimiter,
{{tmax,36,"(Subscript[t, max])"},5,tMax},
{{ymax,75000,"(Subscript[y, max])"},50,100000},

TrackedSymbols:>{pop,logBeta,logMu,logGamma,logDelta,logAlpha,inTau,,iInit,iDelay,tmax,ymax},
SaveDefinitions->True,
ControlPlacement->Top
]
]
closeCellGroup[]

Summary and Future Direction

The method provided can produce a “natural” analytic solution to DDE with a single delay (or advance time) not only for a single nonlinear DDE, but for systems of such DDE. It is expected to be straight forward to expand this approach to systems with multiple DDE and to neutral difference differential equations.

The approach can also work if the natural DDE are expected to minimize the path lengths of the solutions curve, although there are additional challenges with the possible issues of dealing with infinite dimensional metric spaces.

There are certainly possibilities of developing faster and more robust methods for approximating the solution, especially over long ranges of time both forward and backward from the starting times.

Finally, I believe that utilizing this approach to a wide range of existing problems will further show the value of these “natural” analytic solutions to DDE.

Notes on advantages of using Wolfram for this project

- Already has a DDE solver in NDSolve and ParametricNDSolve for specified history
- Has wide range of interpolation functions
- Has wide range of special functions
- Can handle the symbolic computations needed for analysis of different approaches
- Can solve a wide range of problems including implicit functions
- Provides interactive graphic for studying the impact of various inputs
- Provides access to wide range of data sets

References

[AV] J. Arino, P. van den Driessche “Time delays in Epidemic Models; Modeling and Numerical Considerations” in “Delay Differential Equations and Applications”, O. Arino (ed.) Springer, 2006.

[AU] F.M.Asi, A.G.Ulsoy “Analytical Solution of a System of Homogeneous DDE via the Lambert Function”, Proceeding of the American Control Conference 2000.  Author’s Researchgate Link

[BC] Bellman, R., & Cooke, K. L. “Differential-difference equations”. Santa Monica, CA: RAND Corporation, 1963. https://www.rand.org/pubs/reports/R374.html.

[BK] Balachandran, B., Kalmár-Nagy, T., & Gilsinn, D. E. “Delay differential equations”. Recent advances and new directions. Berlin: Springer, 2009. https://link.springer.com/content/pdf/10.1007/978-0-387-85595-0.pdf

[FC] C. Falbo “Analytic and Numeric Solutions to the Delay Differential Equationsy'(t) ==α y(t-δ) ”, Joint Meeting of Northern and Southern California Sections of the MAA. 1995. Falbo’s Sonoma Link

[GK] Gyori, I., & Krisztin, T. “Oscillation results for linear autonomous partial delay differential equations”. Journal of mathematical analysis and applications, 1993. 174(1), 204-217. https://www.sciencedirect.com/science/article/pii/S0022247X8371111X

[HJ] Hale, J. K. “History of delay equations.” In Delay Differential Equations and Applications, pp. 1-28. Springer, Dordrecht, 2006. https://link.springer.com/chapter/10.1007/1-4020-3647-7_1

[KV] V. Kaurov “[UPDATES] Resources For Novel Coronavirus COVID-19”,  https://community.wolfram.com/groups/-/m/t/1872608

[NR1] R. B. Nachbar “Epidemiological Models for Influenza and COVID-19—part 1”, Wolfram, 2020. https://www.wolframcloud.com/obj/rnachbar/Published/EpidemiologicalModelsForInfluenzaAndCOVID-19--part_1.nb

[NR4] R. B. Nachbar “Epidemiological Models for Influenza and COVID-19—part 4”, Wolfram, 2020. https://www.wolframcloud.com/obj/rnachbar/Published/EpidemiologicalModelsForInfluenzaAndCOVID-19--part_4.nb

[PX] Pennec, X. “Probabilities and statistics on riemannian manifolds: A geometric approach” (Doctoral dissertation, RR- 5093, INRIA-00071490). 2004. https://hal.inria.fr/inria-00071490/

[PYZ] L. Peng, W. Yang, D. Zhang, C. Zhuge, L. Hong “Epidemic analysis of COVID-19 in China by dynamical modeling”, Updated June 2020, https://arxiv.org/pdf/2002.06563.pdf

[RR] Rand, R. “Differential-delay equations”. In Luo, Sun (Eds.). Complex Systems (pp. 83-117). Springer, Berlin, Heidelberg, 2011. https://link.springer.com/content/pdf/10.1007%2F978-3-642-17593-0_3.pdf

[RM] Roussel, M. R. “Perturbative-Iterative Computation of Inertial Manifolds of Systems of Delay-Differential Equations with Small Delays”. Algorithms, 2020. 13(9), 209. https://www.mdpi.com/1999-4893/13/9/209

[TL] Teodoro, M. F., Lima, P. M., Ford, N. J., & Lumb, P. M. “Numerical approximation of forward-backward differential equations by a finite element method”. Proceedings of CMMSE, 2009. 1010-1019. https://www.math.tecnico.ulisboa.pt/~plima/papers/Proceedings/3.2.26.pdf

[YS] Yi, Sun. Time-delay systems: Analysis and control using the Lambert W function. Diss. University of Michigan, 2009. https://deepblue.lib.umich.edu/bitstream/handle/2027.42/64756/syjo_1.pdf?sequence=1

[ZCW] Y. Zhou, Z. Chen, X. Wu, Z. Tian, L. Cheng, L. Ye “The Outbreak Evaluation of COVID-19 in Wuhan District of China”, 2020. https://arxiv.org/pdf/2002.09640.pdf

Data Sets

saved data

covidHubeiList={{(DateObject["Wed 22 Jan 2020", DateObject[{2020,1,22},"Day","Gregorian",-4.`]]),444,399,28,17},{(DateObject["Thu 23 Jan 2020", DateObject[{2020,1,23},"Day","Gregorian",-4.`]]),444,399,28,17},{(DateObject["Fri 24 Jan 2020", DateObject[{2020,1,24},"Day","Gregorian",-4.`]]),549,494,31,24},{(DateObject["Sat 25 Jan 2020", DateObject[{2020,1,25},"Day","Gregorian",-4.`]]),761,689,32,40},{(DateObject["Sun 26 Jan 2020", DateObject[{2020,1,26},"Day","Gregorian",-4.`]]),1058,964,42,52},{(DateObject["Mon 27 Jan 2020", DateObject[{2020,1,27},"Day","Gregorian",-4.`]]),1423,1302,45,76},{(DateObject["Tue 28 Jan 2020", DateObject[{2020,1,28},"Day","Gregorian",-4.`]]),3554,3349,80,125},{(DateObject["Wed 29 Jan 2020", DateObject[{2020,1,29},"Day","Gregorian",-4.`]]),3554,3341,88,125},{(DateObject["Thu 30 Jan 2020", DateObject[{2020,1,30},"Day","Gregorian",-4.`]]),4903,4651,90,162},{(DateObject["Fri 31 Jan 2020", DateObject[{2020,1,31},"Day","Gregorian",-4.`]]),5806,5461,141,204},{(DateObject["Sat 1 Feb 2020", DateObject[{2020,2,1},"Day","Gregorian",-4.`]]),7153,6736,168,249},{(DateObject["Sun 2 Feb 2020", DateObject[{2020,2,2},"Day","Gregorian",-4.`]]),11177,10532,295,350},{(DateObject["Mon 3 Feb 2020", DateObject[{2020,2,3},"Day","Gregorian",-4.`]]),13522,12722,386,414},{(DateObject["Tue 4 Feb 2020", DateObject[{2020,2,4},"Day","Gregorian",-4.`]]),16678,15677,522,479},{(DateObject["Wed 5 Feb 2020", DateObject[{2020,2,5},"Day","Gregorian",-4.`]]),19665,18483,633,549},{(DateObject["Thu 6 Feb 2020", DateObject[{2020,2,6},"Day","Gregorian",-4.`]]),22112,20677,817,618},{(DateObject["Fri 7 Feb 2020", DateObject[{2020,2,7},"Day","Gregorian",-4.`]]),24953,23139,1115,699},{(DateObject["Sat 8 Feb 2020", DateObject[{2020,2,8},"Day","Gregorian",-4.`]]),27100,24881,1439,780},{(DateObject["Sun 9 Feb 2020", DateObject[{2020,2,9},"Day","Gregorian",-4.`]]),29631,26965,1795,871},{(DateObject["Mon 10 Feb 2020", DateObject[{2020,2,10},"Day","Gregorian",-4.`]]),31728,28532,2222,974},{(DateObject["Tue 11 Feb 2020", DateObject[{2020,2,11},"Day","Gregorian",-4.`]]),33366,29659,2639,1068},{(DateObject["Wed 12 Feb 2020", DateObject[{2020,2,12},"Day","Gregorian",-4.`]]),33366,29612,2686,1068},{(DateObject["Thu 13 Feb 2020", DateObject[{2020,2,13},"Day","Gregorian",-4.`]]),48206,43437,3459,1310},{(DateObject["Fri 14 Feb 2020", DateObject[{2020,2,14},"Day","Gregorian",-4.`]]),54406,48175,4774,1457},{(DateObject["Sat 15 Feb 2020", DateObject[{2020,2,15},"Day","Gregorian",-4.`]]),56249,49030,5623,1596},{(DateObject["Sun 16 Feb 2020", DateObject[{2020,2,16},"Day","Gregorian",-4.`]]),58182,49847,6639,1696},{(DateObject["Mon 17 Feb 2020", DateObject[{2020,2,17},"Day","Gregorian",-4.`]]),59989,50338,7862,1789},{(DateObject["Tue 18 Feb 2020", DateObject[{2020,2,18},"Day","Gregorian",-4.`]]),61682,50633,9128,1921},{(DateObject["Wed 19 Feb 2020", DateObject[{2020,2,19},"Day","Gregorian",-4.`]]),62031,49665,10337,2029},{(DateObject["Thu 20 Feb 2020", DateObject[{2020,2,20},"Day","Gregorian",-4.`]]),62442,48510,11788,2144},{(DateObject["Fri 21 Feb 2020", DateObject[{2020,2,21},"Day","Gregorian",-4.`]]),62662,48637,11881,2144},{(DateObject["Sat 22 Feb 2020", DateObject[{2020,2,22},"Day","Gregorian",-4.`]]),64084,46439,15299,2346},{(DateObject["Sun 23 Feb 2020", DateObject[{2020,2,23},"Day","Gregorian",-4.`]]),64084,46395,15343,2346},{(DateObject["Mon 24 Feb 2020", DateObject[{2020,2,24},"Day","Gregorian",-4.`]]),64287,45044,16748,2495},{(DateObject["Tue 25 Feb 2020", DateObject[{2020,2,25},"Day","Gregorian",-4.`]]),64786,43252,18971,2563},{(DateObject["Wed 26 Feb 2020", DateObject[{2020,2,26},"Day","Gregorian",-4.`]]),65187,41603,20969,2615},{(DateObject["Thu 27 Feb 2020", DateObject[{2020,2,27},"Day","Gregorian",-4.`]]),65596,39572,23383,2641},{(DateObject["Fri 28 Feb 2020", DateObject[{2020,2,28},"Day","Gregorian",-4.`]]),65914,36829,26403,2682},{(DateObject["Sat 29 Feb 2020", DateObject[{2020,2,29},"Day","Gregorian",-4.`]]),66337,34617,28993,2727},{(DateObject["Sun 1 Mar 2020", DateObject[{2020,3,1},"Day","Gregorian",-4.`]]),66907,32610,31536,2761},{(DateObject["Mon 2 Mar 2020", DateObject[{2020,3,2},"Day","Gregorian",-4.`]]),67103,30366,33934,2803},{(DateObject["Tue 3 Mar 2020", DateObject[{2020,3,3},"Day","Gregorian",-4.`]]),67217,28174,36208,2835},{(DateObject["Wed 4 Mar 2020", DateObject[{2020,3,4},"Day","Gregorian",-4.`]]),67332,25904,38557,2871},{(DateObject["Thu 5 Mar 2020", DateObject[{2020,3,5},"Day","Gregorian",-4.`]]),67466,23972,40592,2902},{(DateObject["Fri 6 Mar 2020", DateObject[{2020,3,6},"Day","Gregorian",-4.`]]),67592,22628,42033,2931},{(DateObject["Sat 7 Mar 2020", DateObject[{2020,3,7},"Day","Gregorian",-4.`]]),67666,21207,43500,2959},{(DateObject["Sun 8 Mar 2020", DateObject[{2020,3,8},"Day","Gregorian",-4.`]]),67707,19486,45235,2986},{(DateObject["Mon 9 Mar 2020", DateObject[{2020,3,9},"Day","Gregorian",-4.`]]),67743,18247,46488,3008},{(DateObject["Tue 10 Mar 2020", DateObject[{2020,3,10},"Day","Gregorian",-4.`]]),67760,16993,47743,3024},{(DateObject["Wed 11 Mar 2020", DateObject[{2020,3,11},"Day","Gregorian",-4.`]]),67773,15593,49134,3046},{(DateObject["Thu 12 Mar 2020", DateObject[{2020,3,12},"Day","Gregorian",-4.`]]),67781,14407,50318,3056},{(DateObject["Fri 13 Mar 2020", DateObject[{2020,3,13},"Day","Gregorian",-4.`]]),67786,13171,51553,3062},{(DateObject["Sat 14 Mar 2020", DateObject[{2020,3,14},"Day","Gregorian",-4.`]]),67790,11755,52960,3075},{(DateObject["Sun 15 Mar 2020", DateObject[{2020,3,15},"Day","Gregorian",-4.`]]),67794,10421,54288,3085},{(DateObject["Mon 16 Mar 2020", DateObject[{2020,3,16},"Day","Gregorian",-4.`]]),67798,9557,55142,3099},{(DateObject["Tue 17 Mar 2020", DateObject[{2020,3,17},"Day","Gregorian",-4.`]]),67799,8685,56003,3111},{(DateObject["Wed 18 Mar 2020", DateObject[{2020,3,18},"Day","Gregorian",-4.`]]),67800,7751,56927,3122},{(DateObject["Thu 19 Mar 2020", DateObject[{2020,3,19},"Day","Gregorian",-4.`]]),67800,6988,57682,3130},{(DateObject["Fri 20 Mar 2020", DateObject[{2020,3,20},"Day","Gregorian",-4.`]]),67800,6285,58382,3133},{(DateObject["Sat 21 Mar 2020", DateObject[{2020,3,21},"Day","Gregorian",-4.`]]),67800,5715,58946,3139},{(DateObject["Sun 22 Mar 2020", DateObject[{2020,3,22},"Day","Gregorian",-4.`]]),67800,5214,59433,3153},{(DateObject["Mon 23 Mar 2020", DateObject[{2020,3,23},"Day","Gregorian",-4.`]]),67800,4768,59879,3153},{(DateObject["Tue 24 Mar 2020", DateObject[{2020,3,24},"Day","Gregorian",-4.`]]),67801,4317,60324,3160},{(DateObject["Wed 25 Mar 2020", DateObject[{2020,3,25},"Day","Gregorian",-4.`]]),67801,3827,60811,3163},{(DateObject["Thu 26 Mar 2020", DateObject[{2020,3,26},"Day","Gregorian",-4.`]]),67801,3431,61201,3169},{(DateObject["Fri 27 Mar 2020", DateObject[{2020,3,27},"Day","Gregorian",-4.`]]),67801,2895,61732,3174},{(DateObject["Sat 28 Mar 2020", DateObject[{2020,3,28},"Day","Gregorian",-4.`]]),67801,2526,62098,3177},{(DateObject["Sun 29 Mar 2020", DateObject[{2020,3,29},"Day","Gregorian",-4.`]]),67801,2049,62570,3182},{(DateObject["Mon 30 Mar 2020", DateObject[{2020,3,30},"Day","Gregorian",-4.`]]),67801,1726,62889,3186},{(DateObject["Tue 31 Mar 2020", DateObject[{2020,3,31},"Day","Gregorian",-4.`]]),67801,1461,63153,3187},{(DateObject["Wed 1 Apr 2020", DateObject[{2020,4,1},"Day","Gregorian",-4.`]]),67802,1283,63326,3193},{(DateObject["Thu 2 Apr 2020", DateObject[{2020,4,2},"Day","Gregorian",-4.`]]),67802,1132,63471,3199},{(DateObject["Fri 3 Apr 2020", DateObject[{2020,4,3},"Day","Gregorian",-4.`]]),67802,987,63612,3203},{(DateObject["Sat 4 Apr 2020", DateObject[{2020,4,4},"Day","Gregorian",-4.`]]),67803,834,63762,3207},{(DateObject["Sun 5 Apr 2020", DateObject[{2020,4,5},"Day","Gregorian",-4.`]]),67803,648,63945,3210},{(DateObject["Mon 6 Apr 2020", DateObject[{2020,4,6},"Day","Gregorian",-4.`]]),67803,577,64014,3212},{(DateObject["Tue 7 Apr 2020", DateObject[{2020,4,7},"Day","Gregorian",-4.`]]),67803,518,64073,3212},{(DateObject["Wed 8 Apr 2020", DateObject[{2020,4,8},"Day","Gregorian",-4.`]]),67803,448,64142,3213},{(DateObject["Thu 9 Apr 2020", DateObject[{2020,4,9},"Day","Gregorian",-4.`]]),67803,401,64187,3215},{(DateObject["Fri 10 Apr 2020", DateObject[{2020,4,10},"Day","Gregorian",-4.`]]),67803,351,64236,3216},{(DateObject["Sat 11 Apr 2020", DateObject[{2020,4,11},"Day","Gregorian",-4.`]]),67803,320,64264,3219},{(DateObject["Sun 12 Apr 2020", DateObject[{2020,4,12},"Day","Gregorian",-4.`]]),67803,303,64281,3219},{(DateObject["Mon 13 Apr 2020", DateObject[{2020,4,13},"Day","Gregorian",-4.`]]),67803,244,64338,3221},{(DateObject["Tue 14 Apr 2020", DateObject[{2020,4,14},"Day","Gregorian",-4.`]]),67803,219,64363,3221},{(DateObject["Wed 15 Apr 2020", DateObject[{2020,4,15},"Day","Gregorian",-4.`]]),67803,179,64402,3222},{(DateObject["Thu 16 Apr 2020", DateObject[{2020,4,16},"Day","Gregorian",-4.`]]),67803,146,64435,3222},{(DateObject["Fri 17 Apr 2020", DateObject[{2020,4,17},"Day","Gregorian",-4.`]]),68128,129,63487,4512},{(DateObject["Sat 18 Apr 2020", DateObject[{2020,4,18},"Day","Gregorian",-4.`]]),68128,122,63494,4512},{(DateObject["Sun 19 Apr 2020", DateObject[{2020,4,19},"Day","Gregorian",-4.`]]),68128,109,63507,4512},{(DateObject["Mon 20 Apr 2020", DateObject[{2020,4,20},"Day","Gregorian",-4.`]]),68128,105,63511,4512},{(DateObject["Tue 21 Apr 2020", DateObject[{2020,4,21},"Day","Gregorian",-4.`]]),68128,102,63514,4512},{(DateObject["Wed 22 Apr 2020", DateObject[{2020,4,22},"Day","Gregorian",-4.`]]),68128,97,63519,4512},{(DateObject["Thu 23 Apr 2020", DateObject[{2020,4,23},"Day","Gregorian",-4.`]]),68128,47,63569,4512},{(DateObject["Fri 24 Apr 2020", DateObject[{2020,4,24},"Day","Gregorian",-4.`]]),68128,23,63593,4512},{(DateObject["Sat 25 Apr 2020", DateObject[{2020,4,25},"Day","Gregorian",-4.`]]),68128,12,63604,4512},{(DateObject["Sun 26 Apr 2020", DateObject[{2020,4,26},"Day","Gregorian",-4.`]]),68128,0,63616,4512}};

covidHubeiListex=Map[Append[#,#[[4]]+#[[5]]]&,covidHubeiList];

covidHubeiList//TableForm[#,TableHeadings->{None,{"Date","Confirmed","Quarantined","Recovered","Deaths"}}]&

Notes:
Quarantine = confirmed - recovered - death
(DateObject["Thu 23 Jan 2020", DateObject[{2020,1,23},"Day","Gregorian",-4.`]]) Wuhan in lockdown
 (DateObject["Sat 15 Feb 2020", DateObject[{2020,2,15},"Day","Gregorian",-4.`]]) changed guidelines for determining COVID to include non DNA
(DateObject["Fri 17 Apr 2020", DateObject[{2020,4,17},"Day","Gregorian",-4.`]]) Additional deaths discovered and reported

HKU data (probably quarantined cases)

covidHubeiEarly=MapIndexed[{DatePlus[DateObject[{2019,12,30}] ,First[#2]],#1}&,{29,35,40,46,53,59,59,59,59,52,44,39,39,41,37,41,48,84,153,226,315,401,466,530,587,619}];

Modules

General

debugPrint

debugPrint checks the debug level and determines if level is high enough to print output

ClearAll[debugPrint]

Options[debugPrint]= {debugOrder-> {None,False,"Fatal","Min",True,Full, All}, deBug-> False, tempPrint-> False};

debugPrint[minDebug_,printInfo_List,Opts:OptionsPattern[]]:= Module[{curDebug,debugLevels,printTemp,x,triggerValue,curValue},
curDebug=OptionValue[deBug];debugLevels=OptionValue[debugOrder];printTemp=OptionValue[ tempPrint];
triggerValue=
(x= Flatten[Position[debugLevels,minDebug]];
If[Length[x]>0,
x[[1]],
Print["Error in debugPrint call, trigger input:",minDebug, ": is not in allowed debug values, has Head:",Head[minDebug]];
0]
);
curValue=
(x= Flatten[Position[debugLevels,curDebug]];
If[Length[x]>0,
x[[1]],
Print["Error in debugPrint call, current debug level:",curDebug,": is not in allowed debug values, has Head:",Head[curDebug]];
0]
);
If[curValue≥ triggerValue,If[printTemp, Apply[PrintTemporary,printInfo],Apply[Print,printInfo]]]
];

closeCellGroup

closeCellGroup[numCells_] :=If[ numCells>=0,(SelectionMove[EvaluationCell[], Next,Cell,numCells]; 
FrontEndExecute[FrontEndToken["SelectionCloseUnselectedCells"]])]

closeCellGroup[] := (SelectionMove[EvaluationCell[], Next,Cell]; 
FrontEndExecute[FrontEndToken["SelectionCloseUnselectedCells"]])

DDE

other Modules

taylorSeries[tsList_List]:=Function[{t},Sum[Evaluate[tsList[[i+2]] ]((((t-Evaluate[tsList[[1]]]))^i)/(i!)),{i,0,Evaluate[Length[tsList]]-2}]];

seriesCoeff[tsSeries_Function,tvar_,tval_,ncoef_]:=Module[{tsValues},Prepend[Table[Flatten[{D[tsSeries[tvar],{tvar,i}]}],{i,0,ncoef}],tvar]/.tvar-> tval]

forwardCalcPtDDE[ddeEq_List,ddePointVal_List,ddeLeadVal_List,vars_List]:=Module[{forCalcPtVal,forCalcOut,tVarName,varNames,tdelay ,delayVars,tdelayVarRules,maxDeriv,forEquations,curValRules,preValRules,forCalcEq},
(* Forward calculate the lead point  derivatives for a single delay τ 
Currently assumes ddePointVal are in the same order as ddeLeadVal list and its values are given
the lead derivatives will will be calculated 
;
__Input______________
vars - a list of variables names with the first variable the independent scaler (t) and the rest the functional names in order used to list equations 
ddeEQ - list of equations in the form D[(Subscript[V, i])[t]]=(Subscript[F, i])[t,(Subscript[V, j])[t],...,(Subscript[V, j])[t-τ],...] in order given by vars[[2;;-1]]
ddePointVal - list with time and variable values and derivatives for each functional  
ddeLeadVal - list with time and variable values at corresponding lead points for each functional  
__Output______________
forCalcOut - List corresponding to ddePointVal with 1 more dimension of derivatives as extention of ddeLeadVal at t+τ
*)
tVarName=First[vars];varNames=Rest[vars];
tdelay = Union[Flatten[Map[Table[Cases[#[[2]],(varNames[[i]][tVarName+td_]->td),{0,Infinity}],{i,1,Length[varNames]}]&,ddeEq]]];
If[Length[tdelay]>1,Print["More tham one time delay is used:",tdelay," only one delay is handled"]];
tdelay =-First[tdelay ];(* tdelay is expected to be the amount of delay and will normally be positive *)
delayVars=Map[(Table[D[varNames[[j]][t],{t,i}],{i,0,Length[#]-3},{j,1,Length[varNames]}]/.{t-> #[[1]]+tdelay})&,ddePointVal];
tdelayVarRules=Map[(#[tVarName+td_]->(Subscript[#, τ])[tVarName])&,varNames];
maxDeriv=Max[Map[Length,ddePointVal]]-2;
forEquations=Table[D[varNames[[j]][tVarName],{tVarName,i+1}]==D[( ddeEq[[j,2]]/.tdelayVarRules),{tVarName,i}],{i,0,maxDeriv},{j,1,Length[varNames]}];
forCalcPtVal=Map[Take[#,2]&,ddeLeadVal];
Do[
curValRules= MapThread[Append[Flatten[Table[ {D[(Subscript[varNames[[j]], τ])[tVarName],{tVarName,i}]->  #1[[2+i,j]],D[varNames[[j]][tVarName],{tVarName,i}]->  #2[[2+i,j]]},{j,1,Length[varNames]},{i,0,k-1}]],tVarName-> #2[[1]]]&,{ddePointVal,forCalcPtVal}];
forCalcEq=Map[(Flatten[forEquations[[k]]/.#])&,curValRules];
forCalcPtVal=MapThread[Append[#1,#2[[All,2]]]&,{forCalcPtVal,forCalcEq}];
,{k,Length[forEquations]}];
(*
forCalcOut=MapThread[{#3[[1]]-tdelay,NSolve[#1,#2]}&,{forCalcEq,forCalcVars,ddePointVal}]; *)
{{delayVars,forEquations,forCalcEq},(* forCalcOut},*)forCalcPtVal}]

splineTaylorPoints Module

splineTaylorPoints generates lowest possible order polynomial spline fit between given points with a specified finite number of derivatives at each point so there is the same continuity at the given points and at least a continuity of the lower endpoint for each knot point created between them. The points are specified by a list with time, value and increasing derivatives specified.

Output will be list of information and the spline function specified as a piecewise function.

This version assumes at most a 1 derivate difference between points, points are specified in time order,  and ...

Options[splineTaylorPoints]= {deBug->False};

splineTaylorPoints[tsValues_List,minKnots_Integer,Opts:OptionsPattern[]]:= Module[{dbugOpts,ckPrint,npoints,tVal,sValues,nDeriv,valuesAsList,funCount,numFun,orderGroup,startOrder,stopOrder,aSpline,nSpace,pOrder,tKnots,tKnotRight, nKnots,lowDeriv,curveSpline,splineConstraints,splineParm,knotValues,curveSplineOut,tsKnots,curveSplineFun,maxExtraKnots},
(* splineTaylorPoints calculates a polynomial spline fit between given points with values and derivatives up to i-th order;
between each pair of points, the spline polynomials have order that is 1 greater than the initial point order;
__Input______________;
  tsValues - a list with the first position as the independent scaler (t) and the rest the value and derivatives at the scalar order 
minKnots - minimum number of knots to be used between values (actual knots may be higher if points have at least minKnots derivatives);
__Output______________;
  Knotspacing - minimum distance between two knots
knotValues - list of "time"/ values/ derivatives at each knot 
curveSplineFun - corresponding spline function ;
*)
dbugOpts=FilterRules[{Opts}, Options[debugPrint]];
ckPrint=OptionValue[deBug];
(* Number of points in list *)
npoints=Length[tsValues];
{tVal,sValues,nDeriv}=Transpose[Map[{First[#],Rest[#],Length[#]-2}&,tsValues]];
debugPrint[True,{TableForm[Prepend[Transpose[{tVal,sValues,nDeriv}],{"time","Value Array", "nDeriv"}]]},deBug-> ckPrint];
(* sValues could be in list of values or just values *)
If[Apply[Plus,Map[Count[#,_List]&,sValues]]>0,
valuesAsList=True;
funCount=Union[Flatten[Map[Map[Length,#]&,sValues]]];
numFun=First[funCount];
If[Length[funCount]>1,Print["Error in splineTaylorPoints input: number of functions in input is not consistence, # of Functions found:",funCount]];
,
valuesAsList=False;numFun=1;
sValues=Map[Transpose,Map[List,sValues]];
(* Put sValues in List form *)
];
debugPrint[True,{"ValueAsList:",valuesAsList," numFun:",numFun," sValues:",sValues},deBug-> ckPrint];
debugPrint[True,{"tValues:",tVal," nDeriv:",nDeriv},deBug-> ckPrint];

{maxExtraKnots,tKnots,knotValues,curveSplineOut}=splineTaylorPoints[tVal,sValues,nDeriv,numFun, minKnots,Opts];

(* tKnotRight (right side of segment) *)
tKnotRight=Append[Drop[tKnots,1],tVal[[-1]]];
debugPrint[True,{{"tKnots:",tKnots},{"tKnotRight:",tKnotRight}},deBug-> ckPrint];
(* Check if Functional list is needed *)
If[Not[valuesAsList], (* remove List from function values *)
knotValues=Flatten[Flatten[knotValues,{3}],1];
curveSplineOut=Flatten[curveSplineOut];
debugPrint[True,{"Flatten list-knotValues:",knotValues," curveSplineOut:",curveSplineOut},deBug-> ckPrint];
];
(* Generate time-values array *)(*
tsKnots=Flatten[Table[Prepend[knotValues[[k,j+1]],tKnots[[k,j+1]]],{k,1,npoints-1},{j,0,nKnots[[k]]}],1];*)
tsKnots=MapThread[{#1,#2}&,{tKnots,knotValues}];
(* Calculate curveSplineFun *)
(*
curveSplineFun=Apply[Piecewise,Append[Prepend[Flatten[MapThread[ Function[{s1,s2,s3},MapThread[  {#1,#2<= t< #3}&,{s1,s2,s3},1]][#1,#2,#3]&,{curveSplineOut,tKnots,tKnotRight},1],1],{curveSplineOut[[1,1]]/.t-> tVal[[1]],t<tVal[[1]]}],{curveSplineOut[[-1,-1]]/. t-> tVal[[-1]],True}]];*)
If[valuesAsList,
curveSplineFun=Map[Piecewise[Append[Prepend[MapThread[ {#1,#2<= t< #3}&,{#,tKnots,tKnotRight}],{#[[1]]/.t-> tVal[[1]],t<tVal[[1]]}],{#[[-1]]/. t-> tVal[[-1]],True}]]&,Transpose[curveSplineOut]]
,
curveSplineFun=Piecewise[Append[Prepend[MapThread[ {#1,#2<= t< #3}&,{curveSplineOut,tKnots,tKnotRight}],{curveSplineOut[[1]]/.t-> tVal[[1]],t<tVal[[1]]}],{curveSplineOut[[-1]]/. t-> tVal[[-1]],True}]]
];
{maxExtraKnots,tsKnots,curveSplineFun}]

splineTaylorPoints[tVal_List,sValues_List,nDeriv_List,numFun_Integer, minKnots_Integer,Opts:OptionsPattern[]]:= Module[{dbugOpts,ckPrint,npoints,dupPosition,dupPos,setCount,setStart,setLen,funCount,orderGroup,startOrder,stopOrder,aSpline,nSpace,nLeftDeriv,nRightDeriv,pOrder,tKnots,tKnotRight,tDistance, nKnots,lowDeriv,curveSpline,splineConstraints,sParm,splineParm,kValues,knotValues,cSplineOut,curveSplineOut,curveSplineFun},
(* splineTaylorPoints calculates a polynomial spline fit between given points with values and derivatives up to i-th order;
between each pair of points, the spline polynomials have order that is 1 greater than the initial point order;
__Input______________;
  tsValues - a list with the first position as the independent scaler (t) and the rest the value and derivatives at the scalar order 
minKnots - minimum number of knots to be used between values (actual knots may be higher if points have at least minKnots derivatives);
__Output______________;
  Knotspacing - minimum distance between two knots
knotValues - list of "time"/ values/ derivatives at each knot 
curveSplineFun - corresponding spline function ;
*)
(* sValues is a list of length npoints
 each element i is a list of numFun functional and each functional at point i has nDeriv[i] derivatives
 *)
dbugOpts=FilterRules[{Opts}, Options[debugPrint]];
ckPrint=OptionValue[deBug];
(* Number of points in independent variable (time) list 
if tVal has duplicates, then need to process in sets without duplicates so tVal is in strictly increasing order in each set *)
npoints=Length[tVal];
tDistance=Flatten[Apply[{Drop[#,1]-Drop[#,-1]}&,{tVal}]];

dupPosition=Flatten[Position[tDistance,n_/;n==0.0]];dupPos=Position[tDistance,n_/;n==0.0];
setCount=Length[dupPosition]+1;
setStart=Prepend[dupPosition+1,1];
setLen=Append[dupPosition,npoints]-(Prepend[dupPosition,0]+1);
debugPrint[True,{"setCount:",setCount," dupPosition:",dupPosition},deBug-> ckPrint];
debugPrint[True,{"setStart:",setStart," setLen:",setLen},deBug-> ckPrint];
nKnots=MapThread[Min[#1,#2]+Abs[Sign[#1-#2]]&,Apply[{Drop[#,1],Drop[#,-1]}&,{nDeriv}]];
nSpace=Map[Max[#,minKnots]+1&,nKnots];
(* there are no knots or spaces between duplicates *)
If[setCount>1,
nKnots=ReplacePart[nKnots,Map[#->0&,dupPosition]];
nSpace=ReplacePart[nSpace,Map[#->0&,dupPosition]];
];
lowDeriv=MapThread[Min[#1,#2]&,Apply[{Drop[#,1],Drop[#,-1]}&,{nDeriv}]];
pOrder=lowDeriv+1;
nLeftDeriv=Drop[nDeriv,-1];
nRightDeriv=Drop[nDeriv,1];
debugPrint[True,{TableForm[Prepend[Transpose[{nLeftDeriv,nRightDeriv,lowDeriv,nKnots,pOrder,nSpace}],{"Left Deriv order","Right Deriv","Low deriv", "nKnots", "Polyorder","nSpace"}]]},deBug-> ckPrint];
(* determine positions of knots (left side of segment) *)
tKnots=Table[tVal[[i]]+If[nSpace[[i]]==0,0,If[j<=nKnots[[i]]/2,j*(tVal[[i+1]]-tVal[[i]])/nSpace[[i]],(nSpace[[i]]-1-nKnots[[i]]+j)*(tVal[[i+1]]-tVal[[i]])/nSpace[[i]]]],{i,1,npoints-1},{j,0,nKnots[[i]]}];
(* tKnotRight (right side of segment) *)
tKnotRight=ReplacePart[tKnots,Map[Rule[#[[1]],#[[2]]]&,Transpose[ReplacePart[MapAt[RotateLeft,Transpose[Flatten[MapIndexed[List[#2,#1]&,tKnots,{2}],1]],-1],{-1,-1}->tVal[[-1]]]]]];
debugPrint[True,{{"tKnots:",tKnots},{"tKnotRight:",tKnotRight},{"tDistance:",tDistance}},deBug-> ckPrint];
(* Create list for output *)
knotValues={};
curveSplineOut={};
(* Splines are split when time is duplicated *) 
Do[
(* create curveSpline functions 
  the spline for each functionals is independent, so can solve curveSpline for each functional separately *)
curveSpline=Table[Sum[(Subscript[aSpline, k,j,i])((((t - tKnots[[k+setStart[[m]]-1,j+1]]))^i)/(i!)),{i,0,pOrder[[k+setStart[[m]]-1]]}],{k,1,setLen[[m]]},{j,0,nKnots[[k+setStart[[m]]-1]]}];
debugPrint[True,{"curveSpline:",curveSpline},deBug-> ckPrint];

splineParm={};
(* Solve for each functional *)
debugPrint[True,{"Set #:",m," Length of set:",setLen[[m]]," Set Start:",setStart[[m]]," order for set:",pOrder[[1+setStart[[m]]-1;;setLen[[m]]+setStart[[m]]-1]]},deBug-> ckPrint];
Do[(* Determine constraints needed for curveSpline for each function *)
splineConstraints=Table[ (* Over each spline piece in current set *)
Join[ (* Left hand side of spline *)
 Table[(D[curveSpline[[k,1]],{t,i}]/.{t->tVal[[k+setStart[[m]]-1]]})==sValues[[k+setStart[[m]]-1,i+1,l]],{i,0,nLeftDeriv[[k+setStart[[m]]-1]]}],
(* Right hand side of spline *)
Table[If[Length[sValues[[k+setStart[[m]]]]]<i+1,debugPrint[True,{"sValues does not exist for point:",tVal[[k+setStart[[m]]]]," deriv:",i+1},deBug-> ckPrint]; Nothing,
(D[curveSpline[[k,-1]],{t,i}]/.{t->tVal[[k+setStart[[m]]]]})==sValues[[k+setStart[[m]],i+1,l]]],{i,0,nRightDeriv[[k+setStart[[m]]-1]]}],
(* Constraints for knots between given times *)Flatten[Table[(D[curveSpline[[k,j]],{t,i}]/.{t->tKnots[[k+setStart[[m]]-1,j+1]]})==(D[curveSpline[[k,j+1]],{t,i}]/.{t->tKnots[[k+setStart[[m]]-1,j+1]]}),{j,1,nKnots[[k+setStart[[m]]-1]]},{i,0,pOrder[[k+setStart[[m]]-1]]-1}]]
],{k,1,setLen[[m]]}];
debugPrint[True,{"Fun:",l," splineConstraints:",splineConstraints},deBug-> ckPrint];
(* Solving for variables in curveSpline to satisfy spline *)
sParm=MapThread[Solve[#1,Cases[Variables[#2],Except[t]]]&,{splineConstraints,curveSpline}];
debugPrint[True,{"sParm:",sParm},deBug-> ckPrint];
splineParm=Append[splineParm,Map[First,sParm]];
debugPrint[True,{"splineParm:",splineParm},deBug-> ckPrint];
,{l,1,numFun}];
(* knotValues groups the function value with each derivative *)
kValues=Table[Table[(Subscript[aSpline, k,j,i])/.splineParm[[l,k]],{i,0,pOrder[[k+setStart[[m]]-1]]},{l,1,numFun}],{k,1,setLen[[m]]},{j,0,nKnots[[k+setStart[[m]]-1]]}];
knotValues=Join[knotValues,kValues];
cSplineOut=Map[Transpose,Transpose[Table[MapThread[(#1/.#2)&,{curveSpline,splineParm[[l]]}],{l,1,numFun}]]];
curveSplineOut=Join[curveSplineOut,cSplineOut];
,{m,1,setCount}];
(* remove duplicate from tKnots and flatten *)
If[setCount>1,
tKnots=Delete[tKnots,dupPos];tKnotRight=Delete[tKnotRight,dupPos];
];
tKnots=Flatten[tKnots];tKnotRight=Flatten[tKnotRight];
debugPrint[True,{{"tKnots:",tKnots},{"tKnotRight:",tKnotRight}},deBug-> ckPrint];
(* Flatten knotValues and curveSplineOut *)
knotValues=Flatten[knotValues,1];curveSplineOut=Flatten[curveSplineOut,1];
debugPrint[True,{"knotValues:",knotValues},deBug-> ckPrint];
debugPrint[True,{"curveSplineOut:",curveSplineOut},deBug-> ckPrint];
(* need to create a list of curveSplineFun for multiple functionals *)
(* To delete:
curveSplineFun=Piecewise[Append[Prepend[Flatten[MapThread[ Function[{s1,s2,s3},MapThread[  {#1,#2<= t< #3}&,{s1,s2,s3},1]][#1,#2,#3]&,{curveSplineOut,tKnots,tKnotRight},1],1],{curveSplineOut[[1,1]]/.t-> tVal[[1]],t<tVal[[1]]}],{curveSplineOut[[-1,-1]]/. t-> tVal[[-1]],True}]];
*)
{Max[nSpace-1],tKnots,knotValues,curveSplineOut}
]

forwardIntegralDDE

ClearAll[forwardIntegralDDE]

Options[forwardIntegralDDE]= {deBug->False,numericCalc->True};

forwardIntegralDDE[ddeEq_List,lagFunction_List,lagRange_List,vars_List,initCond_List,Opts:OptionsPattern[]]:=Module[{dbugOpts,ckPrint,numeric,forCalcPtVal,forCalcOut,tVarName,varNames,numFun,tdelay ,delayFun,lagFunPosition,lagFunReplaceRules,lagDdeEq,forwardEquations},
(* Forward calculate the lead point  derivatives for a single delay τ 
Currently assumes ddePointVal are in the same order as ddeLeadVal list and its values are given
the lead derivatives will will be calculated 
;
__Input______________
vars - a list of variables names with the first variable the independent scaler (t) and the rest the functional names in order used to list equations 
ddeEQ - list of equations in the form D[(Subscript[V, i])[t]]=(Subscript[F, i])[t,(Subscript[V, j])[t],...,(Subscript[V, j])[t-τ],...] in order given by vars[[2;;-1]]
ddePointVal - list with time and variable values and derivatives for each functional  
ddeLeadVal - list with time and variable values at corresponding lead points for each functional  
__Output______________
forCalcOut - List corresponding to ddePointVal with 1 more dimension of derivatives as extention of ddeLeadVal at t+τ
*)
dbugOpts=FilterRules[{Opts}, Options[debugPrint]];
ckPrint=OptionValue[deBug];
numeric=OptionValue[numericCalc];

tVarName=First[vars];varNames=Rest[vars];numFun=Length[varNames];
tdelay = Union[Flatten[Map[Table[Cases[#[[2]],(varNames[[i]][tVarName+td_]->td),{0,Infinity}],{i,1,Length[varNames]}]&,ddeEq]]];
If[Length[tdelay]>1,Print["More tham one time delay is used:",tdelay," only one delay is handled"]];
tdelay =-First[tdelay ];(* tdelay is expected to be the amount of delay and will normally be positive *)

(* Replacing the time delay functional in the DDE with the lagFunction and integrating over lagRange + tdelay *)
delayFun=lagFunction/.{tVarName-> tVarName-tdelay};
lagFunPosition=Map[Table[Position[#,varNames[[i]][tVarName+td_],{0,Infinity}],{i,1,Length[varNames]}]&,ddeEq];
debugPrint[True,{"Lag Function Position:",lagFunPosition},deBug-> ckPrint];
lagFunReplaceRules=Map[Flatten[MapIndexed[#1->delayFun[[First[#2]]]&,#,{2}],1]&,lagFunPosition];
debugPrint[True,{"Replacement rules:",lagFunReplaceRules},deBug-> ckPrint];
lagDdeEq=MapThread[ReplacePart[#1,#2]&,{ddeEq,lagFunReplaceRules}];
debugPrint[True,{"lag DDE:",lagDdeEq},deBug-> ckPrint];
If[numeric,
forwardEquations=NDSolve[Join[lagDdeEq,initCond],varNames,{tVarName,lagRange[[1]]+tdelay,lagRange[[2]]+tdelay}],
forwardEquations=DSolve[Join[lagDdeEq,initCond],varNames,{tVarName,lagRange[[1]]+tdelay,lagRange[[2]]+tdelay}]
];

forwardEquations]

extendKnotsDDE Module

Options[extendKnotsDDE]= {deBug->False,maxPass->100};

extendKnotsDDE[ddeEq_List,currentLinkValues_List,newLocations_List,vars_List,Opts:OptionsPattern[]]:=Module[{mxPass,dbugOpts,ckPrint,newLinkValues,tVarName,varNames,nFun,tdelay ,myLocations,listOfSets,numSets,curTimeList,leadPosition,initPosition,
tVarLead,tVarCur,linkCur,nearCurtVal,nearCurPos,tVarLag,linkLag,nearLeadPos,nearLeadtVal,stepsize,extValList,iLink,tLead,myLeadVal,iPass,extendCurVal,extVar,debug,extendLeadVal,extNearCur,extNearLead,extendNearLeadVal,extCond,extVal,extRules ,curLeadValues,curInitValues,
newLeadValues,newInitValues,newLeadValuesAll,newInitValuesAll,newLeadValuesFull,newInitValuesFull},
(* Given the current linked values calculate values at newlocations that solve DDE
for DDE with a single delay τ 
Currently assumes that the last array in currentLinkValues is the highest linked (lead) value with a current and lag values also in the list
	any additional currentLinkValues will have its matching pair (below current and below lag)
newLocations are in the same order as ddeLeadVal list and will be a list of new times also occurring in pairs (below current and below lag)
new additional derivatives will be calculated for each currentLinkValues along with Taylor series approximations at each new Locations that satisfy the DDE 
;
__Input______________
vars - a list of variables names with the first variable the independent scaler (t) and the rest the functional names in order used to list equations 
ddeEQ - list of equations in the form D[(Subscript[V, i])[t]]=(Subscript[F, i])[t,(Subscript[V, j])[t],...,(Subscript[V, j])[t-τ],...] in order given by vars[[2;;-1]]
currentLinkValues - list with time and variable values and derivatives for each functional that include linked sets 
newLocations - list of time expecting linked times are not included - 
__Output______________
forCalcOut - List corresponding to ddePointVal with 1 more dimension of derivatives as extention of ddeLeadVal at t+τ
*)
dbugOpts=FilterRules[{Opts}, Options[debugPrint]];
ckPrint=OptionValue[deBug];
mxPass=OptionValue[maxPass];
(* Get variables and time delay *)
tVarName=First[vars];varNames=Rest[vars];nFun=Length[varNames];
tdelay = Union[Flatten[Map[Table[Cases[#[[2]],(varNames[[i]][tVarName+td_]->td),{0,Infinity}],{i,1,Length[varNames]}]&,ddeEq]]];
If[Length[tdelay]>1,Print["More than one time delay is used:",tdelay," only one delay is handled"]];
tdelay =-First[tdelay ];(* tdelay is expected to be the amount of delay and will normally be positive *)
(* Determine number of sets of location list are given *)
If[Count[newLocations,_List]==0,
listOfSets=False;myLocations= {newLocations},
listOfSets=True;myLocations= newLocations];
numSets=Length[myLocations];
curTimeList=Map[First,currentLinkValues];
leadPosition=Map[If[#[[1]]>#[[2]],First[Position[curTimeList,First[#]]],Last[Position[curTimeList,First[#]]]]&,myLocations];
initPosition=Map[If[#[[1]]>#[[2]],First[Position[curTimeList,First[#]-tdelay]],Last[Position[curTimeList,First[#]-tdelay]]]&,myLocations];
debugPrint[True,{"numSets:",numSets," curTimeList:",curTimeList," leadPosition:",leadPosition," initPosition:",initPosition},deBug-> ckPrint];

(* Completed variable setup - Above defined variable should be unchanged *)
newInitValuesAll={};newLeadValuesAll={};newInitValuesFull={};newLeadValuesFull={};
Do[ (* over iSet *)
(* Determine initial and lead values for current set of points *)
curLeadValues=currentLinkValues[[First[leadPosition[[iSet]]]]];
curInitValues=currentLinkValues[[First[initPosition[[iSet]]]]];
debugPrint[True,{"curLeadValues:",curLeadValues," curInitValues:",curInitValues},deBug-> ckPrint];
Do[ (* over points in set *)
tVarLead=First[curLeadValues]; 
tVarCur=First[curInitValues];  
(* starting from latest time process newLocations of Set *)
nearLeadtVal=myLocations[[iSet,jPoints]];
nearCurtVal=nearLeadtVal - tdelay;
debugPrint[True,{"tVarLead:",tVarLead," nearLeadtVal:",nearLeadtVal," tVarCur:",tVarCur," nearCurtVal:",nearCurtVal},deBug-> ckPrint];

stepsize=nearLeadtVal-tVarLead;
debugPrint[True,{"stepsize:",stepsize},deBug-> ckPrint];
extValList=extendKnotsDDE[ddeEq,{curInitValues,curLeadValues},stepsize,vars,tdelay,Opts];
debugPrint[True,{"Complete pair call:extValList:",extValList},deBug-> ckPrint];
{extendCurVal,extNearCur,extendLeadVal,extNearLead}=Last[extValList];
debugPrint[True,{"extendLeadVal:",extendLeadVal," extNearLead:",extNearLead},deBug-> ckPrint];

(* build Full list *)
If[stepsize<0,newInitValuesFull=Join[{extNearCur,extendCurVal},newInitValuesFull];newLeadValuesFull=Join[{extNearLead,extendLeadVal},newLeadValuesFull]
,
newInitValuesFull=Join[newInitValuesFull,{extendCurVal,extNearCur}];newLeadValuesFull=Join[newLeadValuesFull,{extendLeadVal,extNearLead}]
];

(* build Set List *)
If[jPoints==2,newInitValues={extendCurVal};newLeadValues={extendLeadVal}];
newInitValues=Append[newInitValues,extNearCur];
newLeadValues=Append[newLeadValues,extNearLead];

(* Change to next point *)
curLeadValues=Drop[extNearLead,-1];
curInitValues=Drop[extNearCur,-1];
debugPrint[True,{"curLeadValues:",curLeadValues,"curInitValues:",curInitValues},deBug-> ckPrint];
,{jPoints,2,Length[myLocations[[iSet]]]}];

If[First[newLeadValues[[1]]]>First[newLeadValues[[2]]],newLeadValues=Reverse[newLeadValues];newInitValues=Reverse[newInitValues]];

newLeadValuesAll=Append[newLeadValuesAll,newLeadValues];
newInitValuesAll=Append[newInitValuesAll,newInitValues];
,{iSet,1,numSets}];
newLinkValues=Map[Drop[#,-1]&,Flatten[Join[newInitValuesAll,newLeadValuesAll],1]];
{{tdelay,newInitValuesAll,newLeadValuesAll}, newLinkValues}]

extendKnotsDDE[ddeEq_List,pairValues_List,stepsize_,vars_List,tdelay_,Opts:OptionsPattern[]]:=Module[{dbugOpts,ckPrint,sortOutput,newLinkValues,tVarName,varNames,nFun,tVarCur,tVarLead,nearCurtVal,nearLeadtVal,
extVar,extendCurVal,otherPossibleVal,extendLeadVal,extNearCur,extNearLead,extendNearLeadVal,extCond,extVal,extRules ,
newLeadValues,newInitValues,newLeadValuesAll,newInitValuesAll},
(* Given the current paired values calculate values at locations shifted by stepsize that solve DDE
for DDE with a single delay τ 

newLocations are in the same order as ddeLeadVal list and will be a list of new times also occurring in pairs (below current and below lag)
new additional derivatives will be calculated for each currentLinkValues along with Taylor series approximations at each new Locations that satisfy the DDE 
;
__Input______________
vars - a list of variables names with the first variable the independent scaler (t) and the rest the functional names in order used to list equations 
ddeEQ - list of equations in the form D[(Subscript[V, i])[t]]=(Subscript[F, i])[t,(Subscript[V, j])[t],...,(Subscript[V, j])[t-τ],...] in order given by vars[[2;;-1]]
pairValues - list with time and variable values and derivatives for each functional that include linked sets 
stepsize - stepsize for nearby points
tdelay - the time delay of DDE 
__Output______________
forCalcOut - List corresponding to ddePointVal with 1 more dimension of derivatives as extention of ddeLeadVal at t+τ
*)
dbugOpts=FilterRules[{Opts}, Options[debugPrint]];
ckPrint=OptionValue[deBug];
sortOutput=False;
(* Get variables names *)
tVarName=First[vars];varNames=Rest[vars];nFun=Length[varNames];
(* get and check times *)
tVarCur=First[pairValues[[1]]];
tVarLead=First[pairValues[[2]]];
nearCurtVal=tVarCur+stepsize;
nearLeadtVal=tVarLead+stepsize;
(* check tdelay = tVarLead - tVarCur : if tdelay is negative then the larger time needs to be first in pairValues *)
If[Not[tVarLead-tVarCur==tdelay],debugPrint[True,{ "Note: time delay given:",tdelay, " not same as difference of times given:",tVarLead-tVarCur},deBug-> ckPrint]];

(* extend derivative at current with unknown *)
extVar=Table[Unique[y],nFun];
extendCurVal=Append[pairValues[[1]],extVar];
(* calculate lead point derivatives from extended current values with DDE *)
{debug,extendLeadVal}=forwardCalcPtDDE[ddeEq,{extendCurVal},{pairValues[[2]]},vars];
extendLeadVal=First[extendLeadVal];(* forwardCalcPtDDE is given the value as a list so output is list *) 

(* use Taylor series to calculate nearby newLocations *)
extNearCur=seriesCoeff[taylorSeries[extendCurVal],tVarName,nearCurtVal,Length[extendCurVal]-2];
extNearLead=seriesCoeff[taylorSeries[extendLeadVal],tVarName,nearLeadtVal,Length[extendLeadVal]-2];
(* calculate near lead point derivatives with DDE *)
{debug,extendNearLeadVal}=forwardCalcPtDDE[ddeEq,{extNearCur},{extNearLead},vars];
extendNearLeadVal=First[extendNearLeadVal];

(* determine extended variable(s) from matching near lead derivative calculation 
 the near lead times and values are the same because of the way constructed *)
extCond=MapThread[#1==#2&,{extNearLead[[3]],extendNearLeadVal[[3]]}];
debugPrint[True,{"extCond:",extCond},deBug-> ckPrint];
extVal=Sort[NSolveValues[extCond,extVar],Norm[#1]<Norm[#2]&]; (* need to sort by other considerations *)
extRules =Map[MapThread[#1->#2&,{extVar,#}]&,extVal];
debugPrint[True,{"extVal:",extVal},deBug-> ckPrint];
 
otherPossibleVal=Drop[Map[{extendCurVal,extNearCur,extendLeadVal,extNearLead}/.#&,extRules],1];
(* Update point values with calculated values *)
extendCurVal=extendCurVal/.First[extRules];
extendLeadVal=extendLeadVal/.First[extRules];
extNearCur=extNearCur/.First[extRules];
extNearLead=extNearLead/.First[extRules];
If[sortOutput,
{otherPossibleVal,Sort[{extendCurVal,extNearCur,extendLeadVal,extNearLead},#1[[1]]<=#2[[1]]&]},
{otherPossibleVal,{extendCurVal,extNearCur,extendLeadVal,extNearLead}}]
]