2
$\begingroup$

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]][[1]] // 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[0] = 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[0] = 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[222]:= General::ivar: 1/1000 is not a valid variable.
During evaluation of In[222]:= General::ivar: 171/20 is not a valid variable.
During evaluation of In[222]:= General::ivar: 8.55` is not a valid variable.
During evaluation of In[222]:= 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`\)}\)\)
Improve this question
$\endgroup$
3
  • $\begingroup$ Is this BVP different of that discussed and solved with NDSolve on scicomp.stackexchange.com/questions/42244/… ? $\endgroup$
    – Alex Trounev
    Dec 13, 2022 at 10:10
  • $\begingroup$ @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. $\endgroup$
    – mathemania
    Dec 13, 2022 at 11:17
  • $\begingroup$ You are right, with FDM this problem can be solved as shown in my answer. $\endgroup$
    – Alex Trounev
    Dec 13, 2022 at 11:51

1 Answer 1

Reset to default
1
$\begingroup$

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]][[1]] // 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[1] = 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[0] = Table[(1 - i) alpha, {i, 1/10, 90/100, (90/100 - 10/100)/(n)}];

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

Visualization

ListLinePlot[Z /. sol, PlotRange -> All]

Figure1

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]][[1]] // 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]]

Figure 2

Improve this answer
$\endgroup$
4
  • $\begingroup$ 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? $\endgroup$
    – mathemania
    Dec 13, 2022 at 16:15
  • $\begingroup$ 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[1]}, 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? $\endgroup$
    – mathemania
    Dec 13, 2022 at 17:09
  • 1
    $\begingroup$ 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. $\endgroup$
    – Alex Trounev
    Dec 13, 2022 at 17:17
  • 1
    $\begingroup$ Use xx = Range[a, b, h]; lst = Table[{xx[[i]], z[i] /. sol}, {i, n + 2}]; zsol = Interpolation[lst]; $\endgroup$
    – Alex Trounev
    Dec 13, 2022 at 17:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.