# How to increase the speed of Collocation method to solve System of differential equations

all. Kindly suggest me a ways to improve the speed of computations.
I recently came to know about collocation methods and their application to find solution for system of differential equations. Here I will describe the procedure in short. Bellow are the system of equations for 2D system. we will assume that we can represent the solution of above differential equations in terms of polynomials of some order n. We will divide the time(t) space into (n-1) so that we can have n point of discrete time, t(i) where i runs from 1 to n. Now we will assume some values of states X1 and X2 at these points. We can construct a polynomial using these data. In this problem I will be using Lagrange interpolation to approximate the states X1 and X2. These approximate states will be represented with X1a,X2a. If we take first derivative for these approximate states (X1a,X2a) and equate it to original state derivatives at the discretized time points.  Since we initially assumed some values of X1,X2 at ti's . Now we will use NMinimize to find the values of X1,X2's at ti's. such that error is minimized. We can simplify this process and represent the final step as Where D_lag is called as D matrix (Lagrange based D Matrix) and X1i is a column matrix of X1-initial guess values at discrete time ti's and X2i is a column matrix of X1-initial guess values at discrete time ti's. f1i and f2i are the value of actual functions evaluated at these discrete time points.

Here is my Mathematica work.

D matrix based on Lagrange can be calculate using this function

dMatrixLagrange[xpoints_] := Block[{ithDLagrangePCoeff,
dpLagrangList, x22, x, n}, n = Length[xpoints];
ithDLagrangePCoeff[x_, j_] :=
Sum[(1/(xpoints[[j]] - xpoints[[i]]))*
Product[(x22 - xpoints[[m1]])/(xpoints[[j]] -
xpoints[[m1]]), {m1, Complement[Range[n], {i, j}]}],
{i, Complement[Range[n], {j}]}] /. x22 -> x;
dpLagrangList[x_] = (ithDLagrangePCoeff[x, #1] & ) /@
Range[n]; (dpLagrangList[#1] & ) /@ xpoints];


System of equations (We will compare Collocation method with NDSolve )

k = 5.;
tf = 3.;
c[(θ1_)?NumericQ, (θ2_)?NumericQ] := NIntegrate[(Sin[θ1 + t]^2.*Cos[θ2 + t])/E^t, {t, 0., tf}];

sol = NDSolve[{
Derivative[θ1][t] == θ2[t]/E^(Sin[t]/10.) + c[θ1[t], θ2[t]],
Derivative[θ2][t] == ((-k)*Sin[θ1[t]])/E^(Sin[t]/10.),
θ1[0.] == 1.,
θ2[0.] == 1.},
{θ1, θ2},
{t, 0., tf}];


Now discretizing the time points. Where npoints - number of time points needed. tpoints- list of time values, dotX- vector consists for first derivative functions

npoints = 11;
tpoints = Subdivide[0.0, tf, npoints - 1];
dotX[(arg_)?(VectorQ[#1, NumericQ] & )] :=Block[{t, θ1, θ2},
t = arg[];
θ1 = arg[];
θ2 = arg[];
{  θ2/E^(Sin[t]/10.) + c[θ1,θ2],
((-k)*Sin[θ1])/E^(Sin[t]/10.)  }
];


Making discrete points symbolically

θ1aData = Table[ToExpression["θ1a" <> ToString[i]], {i, 1,
Length[tpoints]}];

θ2aData = Table[ToExpression["θ2a" <> ToString[i]], {i, 1,
Length[tpoints]}];

1, Length[tpoints]}];


Making LHS side of equations

lhs1 = dMatrixLagrange[tpoints] . θa1θ2aData;


In order to use NMinimize we need objective function in a scalar form. So we take the Norm of LHS-RHS. Final form can be obtained using this function.

eqs[arg_?(VectorQ[#, NumericQ] &)] := Block[{s, sarg},
sarg = Partition[arg, 2];
Sqrt[Total[#^2 & /@ Flatten[(lhs1 /. Thread[s -> arg]) - Table[dotX[Flatten[{tpoints[[i]],
];


Using NMinimize and plotting the results and comparing with NDSolve result.

var = Flatten[θa1θ2aData];
rand = RandomReal[{-2, 2}, Length[var]];
solnmin = NMinimize[{eqs[var], θ1a1 == 1., θ2a1 == 1.}, var];

{  Plot[Evaluate[{θ1[t], θ1poly[t]} /. sol], {t, 0, tf}, PlotStyle -> Automatic,
PlotLegends -> {"θ1", "θ1a"}, AxesLabel -> {"t", "θ's"},  PlotLabel ->
StringJoin["Approximation m points ", ToString[Length[tpoints]], ", nth order ",
ToString[Length[tpoints] - 1]]],
Plot[Evaluate[{θ2[t], θ2poly[t]} /. sol], {t, 0, tf}, PlotStyle -> Automatic,
PlotLegends -> {"θ2", "θ2a"}, AxesLabel -> {"t", "θ's"},
PlotLabel -> StringJoin["Approximation m points ", ToString[Length[tpoints]], ", nth order ",
ToString[Length[tpoints] - 1]]]
}


This approach in Matlab and in Python took few minutes like about 3 minutes to get the result. The way I have written in Mathetmica code might be making the process very slow. to get this result it took me 30 minutes. Kindly suggest some way to speed up the process. Thanking you all.

• Variable var is not defined. Also code looks very unusual for numerical application. Actually it takes 1 s with using NDSolve. What is the advance of your code? – Alex Trounev May 27 at 21:27
• @Alex Trounev, Sorry for the incomplete information. Kindly have a look at the above problem I filled the missing information (var). Yes NDSolve give results with in few seconds this simple system. If the system is very time consuming (for Example-f1,f2 might have to compute NIntegrate based on every time instant). In such cases using collocation methods gives approximate solution with desired accuracy by increasing the number of collocation points (Points where we want to match the approximate solution with exact solution). assuming NMinimize calls less number of time than NDSolve – Gummala Navneeth May 28 at 8:14
• @Alex Trounev , i'm assuming that the way i have written the code might be cause for very slow result. my friends got the results within few seconds in Matlab and Python. – Gummala Navneeth May 28 at 8:18
• @xzczd, I was in an assumption that A.X==B, should be a general form for LinearSolve, and not A.X == f[X] – Gummala Navneeth May 28 at 10:01

First of all, the NDSolve solution can be further improved:

k = 5.;
tf = 3.;
c[θ1_?NumericQ, θ2_?NumericQ] :=
NIntegrate[(Sin[θ1 + t]^2. Cos[θ2 + t])/E^t, {t, 0., tf},
Method -> {Automatic, SymbolicProcessing -> 0}];

sol = NDSolve[{θ1'[t] == θ2[t]/E^(Sin[t]/10.) + c[θ1[t], θ2[t]],
θ2'[t] == -k Sin[θ1[t]]/E^(Sin[t]/10.),
θ1[0.] == 1., θ2[0.] == 1.}, {θ1, θ2}, {t, 0., tf}]; // AbsoluteTiming
(* {0.0730725, Null} *)

ref = sol[[1, All, -1]] // ListLinePlot


Then, as to the implementation of collocation method, NMinimize is undoubtedly a bad choice. Just use FindRoot with the new-in-12.0 "AffineCovariantNewton" method. I've also turned to Chebyshev–Gauss–Lobatto grid to improve the accuracy of the solution:

CGLGrid[{xl_, xr_}, n_] := 1/2 (xl + xr + (xl - xr) Cos[(π Range[0, n - 1])/(n - 1)]);

npoints = 11;
tpoints = CGLGrid[{0., tf}, npoints];
myd = dMatrixLagrange[tpoints];

SetAttributes[c, Listable];
θ1init = {1.}; θ2init = {1.};
tst = FindRoot[{arg1, arg2} |->
With[{t = tpoints, θ1 = θ1init~Join~arg1, θ2 = θ2init~Join~arg2},
Rest /@ {myd . θ1 - (θ2/E^(Sin[t]/10.) + c[θ1, θ2]),
myd . θ2 - (-k Sin[θ1])/E^(Sin[t]/10.)}],
{#, #} &@{ConstantArray[1., npoints - 1]},
Method -> "AffineCovariantNewton"]; // AbsoluteTiming
(* {0.426819, Null} *)

ListPlot[{tpoints, #}\[Transpose] & /@ Join[{θ1init, θ2init}, tst, 2],
PlotRange -> All]~Show~ref ## Remark

1. I've made use of the hidden syntax of FindRoot, a more basic approach is:

  help[arg1_List, arg2_] :=
With[{θ1 = {1.}~Join~arg1, θ2 = {1.}~Join~arg2},
Rest /@ With[{t = tpoints},
{myd . θ1 - (θ2/E^(Sin[t]/10.) + c[θ1, θ2]),
myd . θ2 - (-k Sin[θ1])/E^(Sin[t]/10.)}]]

tst = {arg1, arg2} /.
FindRoot[help[arg1, arg2], {#, ConstantArray[1., npoints - 1]} & /@ {arg1, arg2}]

2. NDSolveFiniteDifferenceDerivative is a better choice to build myd e.g.:

myd = NDSolveFiniteDifferenceDerivative[1, tpoints,
DifferenceOrder -> 2]["DifferentiationMatrix"]

NDSolveFiniteDifferenceDerivative[1, tpoints,
DifferenceOrder -> "Pseudospectral"]["DifferentiationMatrix"] -
dMatrixLagrange[tpoints] // Abs // Max
(* 2.66454*10^-14 *)

3. |-> is introduced in v12.2, if you're not yet in v12.2, use \[Function] instead.

• @ xzczd, link to dMatrixLagrange- " en.wikipedia.org/wiki/Lagrange_polynomial ". in the "Derivatives" section.(For the first derivative, the coefficients are given by--) – Gummala Navneeth May 28 at 10:44
• @ xzczd, hi ,using RootFind approach , I won't be able to give my initial condition (start points of states) .Theta1==1 and Theta2==1. Would you suggest me some way to do that. – Gummala Navneeth May 28 at 10:47
• @Gummala The i.c. is imposed with {1.}~Join~. – xzczd May 28 at 11:28
• @alex Oh I forgot about Chebyshev–Gauss–Lobatto grid. Check my edit. – xzczd May 28 at 12:00
• @GummalaNavneeth Do you mean something like uniform = Range[1., 21]; NDSolveFiniteDifferenceDerivative[1, uniform, DifferenceOrder -> "Pseudospectral", PeriodicInterpolation -> True]["DifferentiationMatrix"]? – xzczd May 28 at 12:55