# Jacobian General::ivar issue

I have a differential equation labelled by eq0 and I want to use finite-difference method to solve it. First, I have to set up the system of equations eqn[i] through a For function. Next, I have to calculate the Jacobian of the coefficient matrix of the system of equations for Newton's method. The parameter n labels the number of equations (or the grid points).

When I calculate the Jacobian up to n=8 there is no problem, but beyond that, say n=10 I encountered an General::ivar issue. I tried calculating the Jacobian using the derivative D function but I encountered the same issue so I guess maybe there is something wrong with how I wrote the code. Any hints?

ClearAll["Global*"]
Needs["VariationalMethods"]
f = 1 - (z[x]/zh)^(d + 1);
L = Sqrt[1 + (z'[x]^2/f)]/z[x]^d;(*Lagrangian*)
eulageq = EulerEquations[L, z[x], x];(*Euler-Lagrange equation*)
s = Solve[eulageq, z''[x]][] // Simplify;(*2nd order EOM*)
eq0 = z''[x] - s[[1, 2]] /. {d -> 3, zh -> 10};

(*Setting up the nonlinear system of equations*)
a = 0;
b = 1;
n = 10;
h = (b - a)/(n + 1);
alpha = 95/10;
beta = 10^-3;
z = alpha;
z[n + 1] = beta;
For[i = 1, i <= n, i++, eqn[i] = Simplify[Collect[eq0 /. {z''[x] -> ((z[i + 1] - 2 z[i] + z[i - 1])/h^2), z'[x] -> ((z[i + 1] - z[i - 1])/(2 h)), z[x] -> z[i]}, z[i]]]; Print["eqn[", i, "] = ", eqn[i]]]

j = 0;
x = Table[(1 - i) alpha, {i, 1/10, 90/100, (90/100 - 10/100)/(n - 1)}];
xr[j] = MapThread[#1 -> #2 &, {Array[z, Length[Table[i, {i, 1, n}]]], x[j]}];
DFx = ResourceFunction["JacobianMatrix"][Table[eqn[i], {i, 1, n}], Table[z[i], {i, 1, n}]] /. xr[j]//N

General::ivar: 1/1000 is not a valid variable.
General::ivar: 171/20 is not a valid variable.
General::ivar: 8.55 is not a valid variable.
General::stop: Further output of General::ivar will be suppressed during this calculation.
JacobianMatrix [{73.26472605,46.03916874,45.22798136,47.78091888,53.1565758,61.92990427,75.89993511,419.515069,262284.3839,-226.5003625},{8.55,7.705555556,6.861111111,6.016666667,5.172222222,4.327777778,3.483333333,2.638888889,1.794444444,0.95}]

(*Jacobian using D function*)
DFx1=D[Table[eqn[i],{i,1,n}],{Table[z[i],{i,1,n}]}]/.xr[j]//N

During evaluation of In:= General::ivar: 1/1000 is not a valid variable.
During evaluation of In:= General::ivar: 171/20 is not a valid variable.
During evaluation of In:= General::ivar: 8.55 is not a valid variable.
During evaluation of In:= General::stop: Further output of General::ivar will be suppressed during this calculation.
\!$$\*SubscriptBox[\(\[PartialD]$$, $${{8.55, 7.705555555555556, 6.861111111111111, 6.016666666666667, 5.1722222222222225, 4.3277777777777775, 3.4833333333333334, 2.638888888888889, 1.7944444444444445, 0.95}}$$]$${73.26472605314348, 46.039168736484946, 45.22798135742989, 47.78091887685099, 53.156575796563885, 61.92990426583507, 75.89993511371034, 419.51506904343574, 262284.38392592594, \(-226.5003624756579$$}\)\)


• Is this BVP different of that discussed and solved with NDSolve on scicomp.stackexchange.com/questions/42244/… ? Dec 13, 2022 at 10:10
• @AlexTrounev It's the same, however, I'm using a different method (finite-difference) to study this BVP since I'll eventually use the best method to study the harder BVP case similar to 275496. I know you've resolved some other BVP issue with wavelet technique but, it doesn't seem to work in the simple case like this. Using shooting method somehow works for the simpler case but I think finite-difference may have an advantage for the harder case for big calculations. Dec 13, 2022 at 11:17
• You are right, with FDM this problem can be solved as shown in my answer. Dec 13, 2022 at 11:51

With a small correction we can compute solution using FindRoot as follows

ClearAll["Global*"]
Needs["VariationalMethods"]
f = 1 - (z[x]/zh)^(d + 1);
L = Sqrt[1 + (z'[x]^2/f)]/z[x]^d;(*Lagrangian*)
eulageq = EulerEquations[L, z[x], x];(*Euler-Lagrange equation*)
s = Solve[eulageq, z''[x]][] // Simplify;(*2nd order EOM*)
eq0 = z''[x] - s[[1, 2]] /. {d -> 3, zh -> 10};

(*Setting up the nonlinear system of equations*)
a = 0;
b = 1;
n = 250;
h = (b - a)/(n + 1);
alpha = 95/10;
beta = 10^-3; Z = Array[z, {n + 2}];
z = alpha;
z[n + 2] = beta; rul =
Table[{z''[x] -> ((z[i + 1] - 2 z[i] + z[i - 1])/h^2),
z'[x] -> ((z[i + 1] - z[i - 1])/(2 h)), z[x] -> z[i]}, {i, 2,
n + 1}];
eqs = Table[eq0 /. rul[[i]], {i, Length[rul]}];

DFx = ResourceFunction["JacobianMatrix"][eqs,
Table[z[i], {i, 2, n + 1}]];

x = Table[(1 - i) alpha, {i, 1/10, 90/100, (90/100 - 10/100)/(n)}];

sol = FindRoot[eqs, Table[{z[i], x[[i]]}, {i, 2, n + 1}],
Jacobian -> DFx, MaxIterations -> 1000, WorkingPrecision -> 30]


Visualization

ListLinePlot[Z /. sol, PlotRange -> All] Using code from my answer here we can compute solution for inverse function $$x(z)$$ as follows

Clear[f, x, z]

f = 1 - (z/zh)^(d + 1);
L = (Sqrt[1 + (x'[z]^(-2)/f)]/z^d) x'[z];(*Lagrangian*)
eulageq = EulerEquations[L, x[z], z];(*Euler-Lagrange equation*)
s = Solve[eulageq, x''[z]][] // Simplify;
eq0 = x''[z] - s[[1, 2]] == 0;

eq1 = eq0 /. {zh -> 10, d -> 3};

zs = 95/10; z0 = 10^-3; sol2 =
NDSolveValue[{eq1, x[zs] == 0, x[z0] == 1}, x, z]


Finally we compare sol and sol2 in one plot

Show[Plot[sol2[z], {z, 10^-3, zs}, PlotRange -> All,
AxesLabel -> {"z", "x"}],
ListPlot[Table[{z[i] /. sol, xx[[i]]}, {i, n + 2}],
PlotStyle -> Red]] • Nice rewrite of the code, I just have one question. Now, instead of using the Jacobian for the Newton's method, you used the Jacobian in FindRoot. From the documentation, it stated that without the Jacobian FindRoot will do extra evaluations for finite differences. Does this mean that FindRoot is using finite difference method internally without specifyng it that specific method? Also, is there some specific reason why you made a shift in the index from i=0 to i=1 and i=n+1 to i=n+2? Dec 13, 2022 at 16:15
• In addition to your code sol, I would like to construct an interpolating function to be substituted back to the Lagrangian and integrate in the domain [x0,xf] (calculate the action). I wrote solution = Join[{z}, Table[sol[[i, 2]], {i, 1, Length[sol]}], {z[n + 2]}]; zfunc = Interpolation[solution, InterpolationOrder -> 3] but zfunc[x] produces an InterpolatingFunction::dmval issue. Why? Dec 13, 2022 at 17:09
• Actually, array Z should be defined before computation. As usual I use origin 1. If you like to use origin 0, just put Z=Array[z,n+1,0], and then rewrite code. Concerning FindRoot , please, note, that option Jacobian ->J is described in tutorial in section Options. Dec 13, 2022 at 17:17
• Use xx = Range[a, b, h]; lst = Table[{xx[[i]], z[i] /. sol}, {i, n + 2}]; zsol = Interpolation[lst]; Dec 13, 2022 at 17:21