# How to speed up FindRoot for nonlinear algebraic equation?

We trying to solve N times N nonlinear algebraic equations using FindRoot in mathematica. Actually its a coupled partial nonlinear differential equation with two variables, we are using pseudospectral method of expanding the function in terms of Chebyshev polynomial and fourier series in u and x respectively.By doing so we are converting Non linear PDE in to algebric equations. Then we evaluate the equations in points of a N times N grid so that the number of variables are equal to number algebric equations. Then we are using Findroot to solve these algebric eqautions.

Our main problem is FindRoot is taking a huge amount of time to run the program when we set the grid size N large like say N=10. Is there a process to speed up FindRoot or any method to parallellize the process of FindRoot.

Clear[Nu, Nx, L, psi, eqz, a1, psi1, z1, ay1]; L = 1.646; Nx = 3; Nu
= 6; mpsi = 0;
mu = 4; b = 0;
psi[x_, u_] =
Sum[psi1[j, k] ChebyshevT[j, 2 u - 1] Cos[2 Pi k x/L], {j, 0,
Nu - 1}, {k, 0, Nx - 1}];
a0[x_, u_] :=
Sum[a1[j, k] ChebyshevT[j, 2 u - 1] Cos[2 Pi k x/L], {j, 0,
Nu - 1}, {k, 0, Nx - 1}];
z[x_, u_] :=
Sum[z1[j, k] ChebyshevT[j, 2 u - 1] Cos[2 Pi k x/L], {j, 0,
Nu - 1}, {k, 0, Nx - 1}];
ay[x_, u_] :=
Sum[ay1[j, k] ChebyshevT[j, 2 u - 1] Cos[2 Pi k x/L], {j, 0,
Nu - 1}, {k, 0, Nx - 1}];
f1 = 1/Sqrt[2]; f2 = 1/Sqrt[2]; psiinfty = Pi/4; psix =
D[psi[x, u], x]; a0x = D[a0[x, u], x]; ayx = D[ay[x, u], x];
zx = D[z[x, u], x]; psiu = D[psi[x, u], u]; a0u =
D[a0[x, u], u]; ayu = D[ay[x, u], u];
zu = D[z[x, u], u];
(*P=1;
G=1;*)
A =
1 + h u^2 psiu^2 + h zu^2 - u^4 a0u^2 +
h u^4 ayu^2; Ax = - a0x^2/h + ayx^2 + zx^2/u^4 +
psix^2/u^2; Axu = - u^4 a0x^2 ayu^2 - zu^2 a0x^2 -
u^2 psiu^2 a0x^2 - u^4 a0u^2 ( zx^2/u^4 + psix^2/u^2 + ayx^2 ) +
h zu^2 ayx^2 + h u^2 psiu^2 ayx^2 + h psiu^2 zx^2/u^2 +
h ayu^2 zx^2 + h zu^2 psix^2 /u^2 + h u^2 ayu^2 psix^2 -
2 h zu psiu zx psix /u^2 - 2 h ayu ayx ( zu zx + u^2 psiu psix ) +
2 u^2 a0u a0x ( u^2 ayu ayx + zu zx /u^2 + psiu psix );
c = psi[x, u] - (1/4) Sin[4 psi[x, u]] -
psiinfty + (1/4) Sin[4 psiinfty]; h = 1 - u^4;
P = Sqrt[u^(-4)  A + Ax +  Axu ];
(*Q = Sqrt[G] Sqrt[r^4 + A + Axu ] - (1/u^2) f1 f2 D[z[x,u],u] - 2 c \
u^2 (a0u  ayx-ayu a0x);*)

G = Sqrt[(f1^2 + 4 (Cos[psi[x, u]])^4 ) (f1^2 +
4 (Sin[psi[x, u]])^4 )];
eqpsi = D[(((2 h psiu)/u^2 - 2 a0x^2 psiu u^2 + 2 ayx^2 h psiu u^2 +
2 a0u a0x psix u^2 - 2 ayu ayx h psix u^2 - (2 h psix zu zx)/
u^2 + (2 h psiu zx^2)/
u^2) Sqrt[(1/2 + 4 Cos[psi[x, u]]^4) (1/2 +
4 Sin[psi[x, u]]^4)])/(2 u^2 \[Sqrt](ayx^2 - a0x^2/h +
psix^2/u^2 - a0x^2 psiu^2 u^2 + ayx^2 h psiu^2 u^2 +
ayu^2 h psix^2 u^2 - a0x^2 ayu^2 u^4 - a0x^2 zu^2 +
ayx^2 h zu^2 + (h psix^2 zu^2)/u^2 + (
1 + h psiu^2 u^2 - a0u^2 u^4 + ayu^2 h u^4 + h zu^2)/u^4 - (
2 h psiu psix zu zx)/u^2 + ayu^2 h zx^2 + zx^2/u^4 + (
h psiu^2 zx^2)/u^2 - 2 ayu ayx h (psiu psix u^2 + zu zx) +
2 a0u a0x u^2 (psiu psix + ayu ayx u^2 + (zu zx)/u^2) -
a0u^2 u^4 (ayx^2 + psix^2/u^2 + zx^2/u^4))), u] +
D[(((2 psix)/u^2 + 2 a0u a0x psiu u^2 - 2 ayu ayx h psiu u^2 -
2 a0u^2 psix u^2 + 2 ayu^2 h psix u^2 + (2 h psix zu^2)/
u^2 - (2 h psiu zu zx)/
u^2) Sqrt[(1/2 + 4 Cos[psi[x, u]]^4) (1/2 +
4 Sin[psi[x, u]]^4)])/(2 u^2 \[Sqrt](ayx^2 - a0x^2/h +
psix^2/u^2 - a0x^2 psiu^2 u^2 + ayx^2 h psiu^2 u^2 +
ayu^2 h psix^2 u^2 - a0x^2 ayu^2 u^4 - a0x^2 zu^2 +
ayx^2 h zu^2 + (h psix^2 zu^2)/u^2 + (
1 + h psiu^2 u^2 - a0u^2 u^4 + ayu^2 h u^4 + h zu^2)/u^4 - (
2 h psiu psix zu zx)/u^2 + ayu^2 h zx^2 + zx^2/u^4 + (
h psiu^2 zx^2)/u^2 - 2 ayu ayx h (psiu psix u^2 + zu zx) +
2 a0u a0x u^2 (psiu psix + ayu ayx u^2 + (zu zx)/u^2) -
a0u^2 u^4 (ayx^2 + psix^2/u^2 + zx^2/u^4))), x] -
2 (-a0x ayu + a0u ayx) u^2 (1 -
Cos[4 psi[x, u]]) + (\[Sqrt](ayx^2 - a0x^2/h + psix^2/u^2 -
a0x^2 psiu^2 u^2 + ayx^2 h psiu^2 u^2 + ayu^2 h psix^2 u^2 -
a0x^2 ayu^2 u^4 - a0x^2 zu^2 + ayx^2 h zu^2 + (
h psix^2 zu^2)/u^2 + (
1 + h psiu^2 u^2 - a0u^2 u^4 + ayu^2 h u^4 + h zu^2)/u^4 - (
2 h psiu psix zu zx)/u^2 + ayu^2 h zx^2 + zx^2/u^4 + (
h psiu^2 zx^2)/u^2 - 2 ayu ayx h (psiu psix u^2 + zu zx) +
2 a0u a0x u^2 (psiu psix + ayu ayx u^2 + (zu zx)/u^2) -
a0u^2 u^4 (ayx^2 + psix^2/u^2 + zx^2/u^4)) (16 Cos[
psi[x, u]] (1/2 + 4 Cos[psi[x, u]]^4) Sin[psi[x, u]]^3 -
16 Cos[psi[x, u]]^3 Sin[
psi[x, u]] (1/2 + 4 Sin[psi[x, u]]^4)))/(2 u^2 Sqrt[(1/2 +
4 Cos[psi[x, u]]^4) (1/2 + 4 Sin[psi[x, u]]^4)]);
eqaz = D[-(1/(
2 u^2)) + ((-2 a0x^2 zu + 2 ayx^2 h zu + (2 h zu)/u^4 + (
2 h psix^2 zu)/u^2 + 2 a0u a0x zx - 2 ayu ayx h zx - (
2 h psiu psix zx)/
u^2) Sqrt[(1/2 + 4 Cos[psi[x, u]]^4) (1/2 +
4 Sin[psi[x, u]]^4)])/(2 u^2 \[Sqrt](ayx^2 - a0x^2/h +
psix^2/u^2 - a0x^2 psiu^2 u^2 + ayx^2 h psiu^2 u^2 +
ayu^2 h psix^2 u^2 - a0x^2 ayu^2 u^4 - a0x^2 zu^2 +
ayx^2 h zu^2 + (h psix^2 zu^2)/u^2 + (
1 + h psiu^2 u^2 - a0u^2 u^4 + ayu^2 h u^4 + h zu^2)/
u^4 - (2 h psiu psix zu zx)/u^2 + ayu^2 h zx^2 + zx^2/
u^4 + (h psiu^2 zx^2)/u^2 -
2 ayu ayx h (psiu psix u^2 + zu zx) +
2 a0u a0x u^2 (psiu psix + ayu ayx u^2 + (zu zx)/u^2) -
a0u^2 u^4 (ayx^2 + psix^2/u^2 + zx^2/u^4))), u] +
D[((2 a0u a0x zu - 2 ayu ayx h zu - (2 h psiu psix zu)/u^2 -
2 a0u^2 zx + 2 ayu^2 h zx + (2 zx)/u^4 + (2 h psiu^2 zx)/
u^2) Sqrt[(1/2 + 4 Cos[psi[x, u]]^4) (1/2 +
4 Sin[psi[x, u]]^4)])/(2 u^2 \[Sqrt](ayx^2 - a0x^2/h +
psix^2/u^2 - a0x^2 psiu^2 u^2 + ayx^2 h psiu^2 u^2 +
ayu^2 h psix^2 u^2 - a0x^2 ayu^2 u^4 - a0x^2 zu^2 +
ayx^2 h zu^2 + (h psix^2 zu^2)/u^2 + (
1 + h psiu^2 u^2 - a0u^2 u^4 + ayu^2 h u^4 + h zu^2)/u^4 - (
2 h psiu psix zu zx)/u^2 + ayu^2 h zx^2 + zx^2/u^4 + (
h psiu^2 zx^2)/u^2 - 2 ayu ayx h (psiu psix u^2 + zu zx) +
2 a0u a0x u^2 (psiu psix + ayu ayx u^2 + (zu zx)/u^2) -
a0u^2 u^4 (ayx^2 + psix^2/u^2 + zx^2/u^4))), x];

eqay = D[((2 ayu h + 2 ayu h psix^2 u^2 - 2 a0x^2 ayu u^4 +
2 a0u a0x ayx u^4 + 2 ayu h zx^2 -
2 ayx h (psiu psix u^2 + zu zx)) Sqrt[(1/2 +
4 Cos[psi[x, u]]^4) (1/2 +
4 Sin[psi[x, u]]^4)])/(2 u^2 \[Sqrt](ayx^2 - a0x^2/h +
psix^2/u^2 - a0x^2 psiu^2 u^2 + ayx^2 h psiu^2 u^2 +
ayu^2 h psix^2 u^2 - a0x^2 ayu^2 u^4 - a0x^2 zu^2 +
ayx^2 h zu^2 + (h psix^2 zu^2)/u^2 + (
1 + h psiu^2 u^2 - a0u^2 u^4 + ayu^2 h u^4 + h zu^2)/
u^4 - (2 h psiu psix zu zx)/u^2 + ayu^2 h zx^2 + zx^2/
u^4 + (h psiu^2 zx^2)/u^2 -
2 ayu ayx h (psiu psix u^2 + zu zx) +
2 a0u a0x u^2 (psiu psix + ayu ayx u^2 + (zu zx)/u^2) -
a0u^2 u^4 (ayx^2 + psix^2/u^2 + zx^2/u^4))) -
2 a0x u^2 (-psiinfty + psi[x, u] + 1/4 Sin[4 psiinfty] -
1/4 Sin[4 psi[x, u]]), u] +
D[((2 ayx + 2 ayx h psiu^2 u^2 + 2 a0u a0x ayu u^4 -
2 a0u^2 ayx u^4 + 2 ayx h zu^2 -
2 ayu h (psiu psix u^2 + zu zx)) Sqrt[(1/2 +
4 Cos[psi[x, u]]^4) (1/2 +
4 Sin[psi[x, u]]^4)])/(2 u^2 \[Sqrt](ayx^2 - a0x^2/h +
psix^2/u^2 - a0x^2 psiu^2 u^2 + ayx^2 h psiu^2 u^2 +
ayu^2 h psix^2 u^2 - a0x^2 ayu^2 u^4 - a0x^2 zu^2 +
ayx^2 h zu^2 + (h psix^2 zu^2)/u^2 + (
1 + h psiu^2 u^2 - a0u^2 u^4 + ayu^2 h u^4 + h zu^2)/
u^4 - (2 h psiu psix zu zx)/u^2 + ayu^2 h zx^2 + zx^2/
u^4 + (h psiu^2 zx^2)/u^2 -
2 ayu ayx h (psiu psix u^2 + zu zx) +
2 a0u a0x u^2 (psiu psix + ayu ayx u^2 + (zu zx)/u^2) -
a0u^2 u^4 (ayx^2 + psix^2/u^2 + zx^2/u^4))) +
2 a0u u^2 (-psiinfty + psi[x, u] + 1/4 Sin[4 psiinfty] -
1/4 Sin[4 psi[x, u]]), x];

eqa0 = D[((-2 a0u +
2 a0x u^2 (psiu psix + ayu ayx u^2 + (zu zx)/u^2) -
2 a0u u^4 (ayx^2 + psix^2/u^2 + zx^2/u^4)) Sqrt[(1/2 +
4 Cos[psi[x, u]]^4) (1/2 +
4 Sin[psi[x, u]]^4)])/(2 u^2 \[Sqrt](ayx^2 - a0x^2/h +
psix^2/u^2 - a0x^2 psiu^2 u^2 + ayx^2 h psiu^2 u^2 +
ayu^2 h psix^2 u^2 - a0x^2 ayu^2 u^4 - a0x^2 zu^2 +
ayx^2 h zu^2 + (h psix^2 zu^2)/u^2 + (
1 + h psiu^2 u^2 - a0u^2 u^4 + ayu^2 h u^4 + h zu^2)/
u^4 - (2 h psiu psix zu zx)/u^2 + ayu^2 h zx^2 + zx^2/
u^4 + (h psiu^2 zx^2)/u^2 -
2 ayu ayx h (psiu psix u^2 + zu zx) +
2 a0u a0x u^2 (psiu psix + ayu ayx u^2 + (zu zx)/u^2) -
a0u^2 u^4 (ayx^2 + psix^2/u^2 + zx^2/u^4))) +
2 ayx u^2 (-psiinfty + psi[x, u] + 1/4 Sin[4 psiinfty] -
1/4 Sin[4 psi[x, u]]), u] +
D[((-((2 a0x)/h) - 2 a0x psiu^2 u^2 - 2 a0x ayu^2 u^4 -
2 a0x zu^2 +
2 a0u u^2 (psiu psix + ayu ayx u^2 + (zu zx)/u^2)) ((1/2 +
4 Cos[psi[x, u]]^4) (1/2 + 4 Sin[psi[x, u]]^4))^(
1/4))/(2 u^2 \[Sqrt](ayx^2 - a0x^2/h + psix^2/u^2 -
a0x^2 psiu^2 u^2 + ayx^2 h psiu^2 u^2 +
ayu^2 h psix^2 u^2 - a0x^2 ayu^2 u^4 - a0x^2 zu^2 +
ayx^2 h zu^2 + (h psix^2 zu^2)/u^2 + (
1 + h psiu^2 u^2 - a0u^2 u^4 + ayu^2 h u^4 + h zu^2)/
u^4 - (2 h psiu psix zu zx)/u^2 + ayu^2 h zx^2 + zx^2/
u^4 + (h psiu^2 zx^2)/u^2 -
2 ayu ayx h (psiu psix u^2 + zu zx) +
2 a0u a0x u^2 (psiu psix + ayu ayx u^2 + (zu zx)/u^2) -
a0u^2 u^4 (ayx^2 + psix^2/u^2 + zx^2/u^4))) -
2 ayu u^2 (-psiinfty + psi[x, u] + 1/4 Sin[4 psiinfty] -
1/4 Sin[4 psi[x, u]]), x];
bc1 = psi[x, 0] - psiinfty;
bc2 = D[psi[x, u], u] /. u -> 0 - mpsi;
bc3 = a0[x, 0] - mu;
bc4 = a0[x, 1];
bc5 = z[x, 0];
bc6 = ay[x, 0] - b x;

tb = Table[{u -> a, x -> a2}, {a, .1, .5, .2}, {a2, .1, 1, .3}];
tbu = tb[[1, 1]][[1, 1]] ->
Flatten[Table[Table[tb[[j, i]][[1, 2]], {i, 1, 4}], {j, 1, 3}]];
tbx = tb[[1, 1]][[2, 1]] ->
Flatten[Table[Table[tb[[j, i]][[2, 2]], {i, 1, 4}], {j, 1, 3}]];

cpt = {tbu, tbx};
dat1 = Flatten[
Table[{psi1[j, k] -> 0, a1[j, k] -> 0}, {j, 0, Nu - 1}, {k, 0,
Nx - 1}]];
dat2 = Flatten[
Table[{ay1[j, k] -> 0, z1[j, k] -> 0}, {j, 0, Nu - 1}, {k, 0,
Nx - 1}]];
dat3 = Flatten[{dat1, dat2}];
Length[dat3]
eqs = {{eqa0, eqay, eqaz, eqpsi} /.
cpt, {bc1, bc2, bc3, bc4, bc5, bc6} /.
x -> Table[tbx[[2, i]], {i, 1, 4}]};
Length[Flatten[eqs]]

72

72

s1 = FindRoot[Re[eqs], dat3 /. Rule -> List,
Method -> {"Newton", "StepControl" -> "TrustRegion"}] //
AbsoluteTiming
Sort[Flatten[eqs /. s1[[2, All]]]] // AbsoluteTiming

• Do you happen to minimize something here? Some background on the equations might help. You may try to compile the equations and their Jacobian in order to speed up evaluation, and, utilizing the Jacobian option of FindRoot. But your equations are really obfuscated and I doubt that Compile will manage to create efficient code from them... – Henrik Schumacher Apr 14 at 9:41
• Moreover, it would be great to see the PDE. – Henrik Schumacher Apr 14 at 9:42
• @Henrik The PDE is given in the above code as eqpsi, eqay, eqaz and eqa0. – Nishal Rai Apr 14 at 9:45
• As general suggestion: For efficiency, you will be forced to shift large parts of the symbolic code to vectorized numerical code. The key to numerical methods involving Chebyshev polynomials is to evalutate the polynomials on an unisolvent set of points and to use the the vectors of function values as degrees of freedom. – Henrik Schumacher Apr 14 at 9:59
• In general, numeric evaluation in machine precision is orders of magnitudes faster than symbolic evaluations. I am not familiar with the Chebyshev pseudospectral method but I am pretty sure that it is not applied like this. – Henrik Schumacher Apr 14 at 10:29