# How to realize NDSolve by NSolve with the implicit constraint equations?

Daniel post the amazing code to realize numerical integrators, I think it is useful to the kind of one-step explicit equation like x[k+1]-x[k]== constant_number. In the general situation, the integrator may be implicit, eg.

deltat = 1*^-3;
integrator[x_, y_,
k_] := {x[[k + 1]] - x[[k]] == deltat*(x[[k + 1]]*x[[k]] + 1),
y[[k + 1]] - y[[k]] == deltat*(y[[k + 1]]*y[[k]]) + 2};


Integrator describe the discrete trajectories of some ODES, after define the initial condition like x[[1]]=0,y[[1]]=0, we can calculate the next point by Solve the pre-defined function integrator.

x = {0, x1};
y = {0, y1};
Solve[integrator[x, y, 1], {x[[2]], y[[2]]}]


After Enter+Shift, output is {{x1 -> 1/1000, y1 -> 2}}

How to compute the next $$k-1$$ point(Unlike above code only can calculate the second point of trajectory), compute points step by step and substitute in step by step seems not smart.

• Your example is not complicated enough for a real implicit scheme. You may solve this equations in every step to see how an implicit method works. However, the equation in this simple example may be solved once for all because the solution does NOT depend on k and can be solved analytically: {xnew -> (deltat + xold)/(1 - deltat xold), ynew -> (-2 - yold)/(-1 + deltat yold)} Nov 17 '21 at 14:15
• @DanielHuber Right, could I just assume the example is complicated, I think the general method is use Solve to find the root of ‘integrators’, I have success to compute the next point of initial condition, the problem maybe to find a loop to make Solve automatically compute more numerical point of ODEs.
– Ben
Nov 17 '21 at 14:39
• If you writeintegrate so that it takes xold,yold, uses NSolve (because it can not be done analytically by Solve ) to get xnew,ynew and then repeat e.g. in a loop, feeding the new coord. as input to integrate. More elegant, you could use NestList. Nov 17 '21 at 14:47
• Assume x[k+1] can’t explicitly express by x[k], so in the general situation, taking ‘Solve’ is a good choice.
– Ben
Nov 17 '21 at 14:48
• @DanielHuber Let me try it, I haven’t use NestList before.
– Ben
Nov 17 '21 at 14:52

I did something similar with stormerVerlet here: Plotting elliptical orbits using Verlet method.

If we break down the operation,

trajectory = NestList[stepIntegrator, condition, nsteps]


where the input and output of stepIntegrator is a point of the trajectory. The paradigm would be of the form

{x1, y1} = myIntegrator[param1, param2,..][{x0, y0}]


In practice, it is helpful to keep the value of the derivative computed at each step, in which case the paradigm becomes

{x1, y1, yp1} = myIntegrator[param1, param2,..][{x0, y0, yp0}]


In these cases, we have stepIntegrator = myIntegrator[param1, param2,...] where the parameters may be omitted or used to define a general method that may depend on step size, the equation, etc.

Example: Implicit Euler. (It could be made more efficient, I suppose, but it demonstrates the idea.)

ClearAll[myImplEuler, myIntegrator];
myImplEuler[
eqn_,       (* ode in terms of x, y, y' *)
{x_, y_},   (* symbolic variables *)
{x0_?NumericQ, y0_?NumericQ, yp0_?NumericQ}, (* current values *)
h_          (* step size *)
] := {y[0], y[1]} /.
FindRoot[{eqn, y == y0 + h*y'} /. {x -> x0 + h, y' -> y[1],
y -> y[0]}, {{y[0], y0}, {y[1], yp0}}];
myIntegrator[
eqn_,   (* ode in terms of x, y, y' *)
vars_,  (* {x, y} *)
h_      (* step size *)
][{x0_, y0_, yp0_}] := {x0 + h,
Sequence @@ myImplEuler[eqn, vars, {x0 + h, y0, yp0}, h]};

stepIntegrator = myIntegrator[y' == y, {x, y}, 0.1];
NestList[stepIntegrator,
{0., 1., 1.},  (* need to provide initial value for y' at the start *)
10][[All, ;; 2]]
ListLinePlot[
{NDSolveValue[{y'[x] == y[x], y[0] == 1}, y, {x, 0, 1}],
%},
PlotLegends -> {NDSolve, myIntegrator}, Epilog -> {Black, Point[%]}]
(* trajectory:
{{0., 1.}, {0.1, 1.11111}, {0.2, 1.23457}, {0.3, 1.37174}, {0.4,
1.52416}, {0.5, 1.69351}, {0.6, 1.88168}, {0.7, 2.09075}, {0.8,
2.32306}, {0.9, 2.58117}, {1., 2.86797}}
*)


stepIntegrator = myIntegrator[y' == y, {x, y}, 2.^-11];
Nest[stepIntegrator, {0., 1., 1.}, 2^11] // Most

(*  {1., 2.71895}  *)


Note the initial value for y' could be solved for by adding such a step to the code. It's used here only as an initial guess for the value of y' in FindRoot at the next step. An accurate value is important only if the ODE eqn has multiple solutions or is likely to fail to converge if the starting point is far away from the true value.

• Fantastic, the idea is clear, now i know how to do it, thank you!
– Ben
Nov 18 '21 at 2:49
• I have some troubles, I have a new idea, use NestList to do the compute and feed process, it seems more clear, but I fails in the last line code(Se my Update answer).
– Ben
Nov 18 '21 at 16:38

Update(this answer have some bugs, but I find some way to compute the result, I will keep update until success. Thanks to everyone)

I Have do it(See the follows code), thanks to Michael and Syed!

deltat = 1*^-3;

integrator[x_, y_, {x0_, y0_}] := {x[1], y[1]} /.
FindRoot[{x[1] - x[0] == deltat,
y[1] - y[0] == deltat} /. {x[0] -> x0, y[0] -> y0}, {{x[1],
x0}, {y[1], y0}}];
integrator[x, y, {0, 0}]
myintegrator[x_, y_][{x0_, y0_}] := integrator[x, y, {x0, y0}];
stepintegrator = myintegrator[x, y]

myintegrator[x, y][{0, 0}]
NestList[stepintegrator, {0, 0}, 3]