2
$\begingroup$

We are trying to solve a large system of implicit system of differential equations for complex variables $\boldsymbol{\alpha}\left(t\right)$ and $\boldsymbol{\lambda}\left(t\right)$. System consists of two types of equations. The first equation (eq1) is $$\sum_{j}\dot{\alpha}_{j,n}S_{ij}+\sum_{j}\alpha_{j,n}S_{ij}\sum_{rh}\left(\dot{\lambda}_{j,rh}L_{ijrh}-\dot{\lambda}_{j,rh}^{\star}\frac{1}{2}\lambda_{j,rh}\right)=\Gamma_{in}\,$$

for every pair of free $i,n$ indices, with $\dot{x}\left(t\right)$ being a time derivative of $x\left(t\right)$, and the second equation (eq2) is

$$\sum_{j,n}\dot{\alpha}_{j,n}\alpha_{i,n}^{\star}S_{ij}\lambda_{j,kq}+\sum_{j,n}\alpha_{i,n}^{\star}\alpha_{j,n}S_{ij}\left(\dot{\lambda}_{j,kq}+\sum_{rh}\dot{\lambda}_{j,rh}L_{ijrh}\lambda_{j,kq}-\sum_{rh}\dot{\lambda}_{j,rh}^{\star}\frac{1}{2}\lambda_{j,rh}\lambda_{j,kq}\right)=\Theta_{ikq}\ ,$$

for every triple of $i,k,q$ indices.

The typical dimensions of $\boldsymbol{\alpha}\left(t\right)$ is $10\times5$ and $10\times50\times5$ for $\boldsymbol{\lambda}\left(t\right)$, so the total number of equations is in the thousands. The typical initial variable values are $\boldsymbol{\alpha}\left(0\right)=1$ and $\boldsymbol{\lambda}\left(0\right)=0$ for all indices, making the system of equations stiff, however, solving the system remains problematic even for small dimensions:

nEL = 3;
nMUL = 2;
nBAT = nEL;
nMOD = 1;

Defining tables of variables we would like to solve for:

a = Table[\[Alpha][i, n][t], {i, nMUL}, {n, nEL}];
l = Table[\[Lambda][i, r, h][t], {i, nMUL}, {r, nBAT}, {h, nMOD}];
ac = Table[\[Alpha]c[i, n][t], {i, nMUL}, {n, nEL}];
lc = Table[\[Lambda]c[i, r, h][t], {i, nMUL}, {r, nBAT}, {h, nMOD}];

and auxiliary variables and functions needed to construct eq1 and eq2:

\[Gamma] = Table[\[CapitalGamma][i, n], {i, nMUL}, {n, nEL}]; 
\[Theta] = Table[\[CapitalTheta][i, k, q], {i, nMUL}, {k, nBAT}, {q, nMOD}];
\[Omega][r_Integer?Positive, h_Integer?Positive] := 25*h /; (r <= nBAT) && (h <= nMOD)
g[n_Integer?Positive, r_Integer?Positive, h_Integer?Positive] := (1/50)*n /; (n <= nEL) && (r <= nBAT) && (h <= nMOD)
J[n_Integer?Positive, m_Integer?Positive] := 1 + KroneckerDelta[n, m]*(5*n - 1) /; (n <= nEL) && (m <= nEL)
S[i_Integer?Positive, j_Integer?Positive] := Exp[Sum[\[Lambda]c[i, k, q][t]*\[Lambda][j, k, q][t] - (1/2)*(Abs[\[Lambda][i, k, q][t]] + Abs[\[Lambda][j, k, q][t]]), {k, nBAT}, {q, nMOD}]] /; (i <= nMUL) && (j <= nMUL)
L[i_Integer?Positive, j_Integer?Positive, r_Integer?Positive,h_Integer?Positive] := \[Lambda]c[i, r, h][t] - (1/2)*\[Lambda]c[j, r, h][t] /; (i <= nMUL) && (j <= nMUL) && (r <= nBAT) && (h <= nMOD)
A[i_Integer?Positive, j_Integer?Positive, r_Integer?Positive, h_Integer?Positive] := \[Omega][r, h]*\[Lambda]c[i, r, h][t]*\[Lambda][j, r, h][t] /; (i <= nMUL) && (j <= nMUL) && (r <= nBAT) && (h <= nMOD)
B[i_Integer?Positive, j_Integer?Positive, r_Integer?Positive,h_Integer?Positive] := -\[Omega][r,h]*(\[Lambda]c[i, r, h][t] + \[Lambda][j, r, h][t]) /; (i <= nMUL) && (j <= nMUL) && (r <= nBAT) && (h <= nMOD)
\[CapitalGamma][i_Integer?Positive,n_Integer?Positive] := -I*Sum[Sum[\[Alpha][j, m][t]*S[i, j]*J[n, m], {m,nEL}] + \[Alpha][j, n][t]*S[i, j]*Sum[(A[i, j, r, h] + g[n, r, h]*B[i, j, r, h]), {r, nBAT}, {h,nMOD}], {j, nMUL}] /; (i <= nMUL) && (n <= nEL)
\[CapitalTheta][i_Integer?Positive, k_Integer?Positive,q_Integer?Positive] := -I*Sum[\[Alpha]c[i, n][t]*\[Alpha][j, m][t]*S[i, j]*J[n, m]*\[Lambda][j, k, q][t], {j, nMUL}, {n, nEL}, {m, nEL}] - I*Sum[\[Alpha]c[i, n][t]*\[Alpha][j, n][t]*S[i, j]*(\[Omega][k,q]*(\[Lambda][j, k, q][t] - g[n, k, q]) + \[Lambda][j, k,q][t]*Sum[(A[i, j, r, h] + g[n, r, h]*B[i, j, r, h]), {r,nBAT}, {h, nMOD}]), {j, nMUL}, {n, nEL}] /; (i <= nMUL) && (k <= nBAT) && (q <= nMOD) 
\[CapitalGamma]c[i_Integer?Positive, n_Integer?Positive] := I*Sum[Sum[\[Alpha]c[j, m][t]*S[j, i]*J[n, m], {m,nEL}] + \[Alpha]c[j, n][t]*S[j, i]*Sum[(A[j, i, r, h] + g[n, r, h]*B[j, i, r, h]), {r, nBAT}, {h,nMOD}], {j, nMUL}] /; (i <= nMUL) && (n <= nEL)
\[CapitalTheta]c[i_Integer?Positive, k_Integer?Positive,q_Integer?Positive] := I*Sum[\[Alpha][i, n][t]*\[Alpha]c[j, m][t]*S[j, i]*J[n, m]*\[Lambda]c[j, k, q][t], {j, nMUL}, {n, nEL}, {m, nEL}] +I*Sum[\[Alpha][i, n][t]*\[Alpha]c[j, n][t]*S[j, i]*(\[Omega][k,q]*(\[Lambda]c[j, k, q][t] - g[n, k, q]) + \[Lambda]c[j, k, q][t]*Sum[(A[j, i, r, h] + g[n, r, h]*B[j, i, r, h]), {r,nBAT}, {h, nMOD}]), {j, nMUL}, {n, nEL}] /; (i <= nMUL) && (k <= nBAT) && (q <= nMOD)

Constructing tables of eq1 for indices $i,n$ and eq2 for indices $i,k,q$:

eq1 = Flatten[Table[Sum[D[\[Alpha][j, n][t], t]*S[i, j], {j, nMUL}] + Sum[\[Alpha][j, n][t]*S[i, j]*(Sum[D[\[Lambda][j, r, h][t], t]*L[i, j, r, h] - (1/2)*D[\[Lambda]c[j, r, h][t], t]*\[Lambda][j, r, h][t], {r, nBAT}, {h, nMOD}]), {j, nMUL}] == \[CapitalGamma][i, n], {i, nMUL}, {n, nEL}]];
eq2 = Flatten[Table[Sum[D[\[Alpha][j, n][t], t]*\[Alpha]c[i, n][t]*S[i, j]*\[Lambda][j, k, q][t] + \[Alpha]c[i, n][t]*\[Alpha][j, n][t]*S[i, j]*(D[\[Lambda][j, k, q][t], t] + Sum[D[\[Lambda][j, r, h][t], t]*L[i, j, r, h]*\[Lambda][j, k, q][t] - (1/2)*D[\[Lambda]c[j, r, h][t], t]*\[Lambda][j, r, h][t]*\[Lambda][j, k, q][t], {r, nBAT}, {h, nMOD}]), {j, nMUL}, {n, nEL}] == \[CapitalTheta][i, k, q], {i, nMUL}, {k, nBAT}, {q, nMOD}]];

Constructing a system to be solved. In order to have an equal number of equations and unknowns $\alpha(t), \lambda(t), \alpha^*(t), \lambda^*(t)$, we additionally have to take complex conjugate of eq1 and eq2:

allEqs = Join[eq1, eq2, eq1 /. {\[Alpha] -> \[Alpha]c, \[Alpha]c -> \[Alpha], \[Lambda] -> \[Lambda]c, \[Lambda]c -> \[Lambda], \[CapitalGamma] -> \[CapitalGamma]c}, eq2 /. {\[Alpha] -> \[Alpha]c, \[Alpha]c -> \[Alpha], \[Lambda] -> \[Lambda]c, \[Lambda]c -> \[Lambda], \[CapitalTheta] -> \[CapitalTheta]c}];

Defining a list of initial values. Again, doubling the number of values by taking a complex conjugate:

initCond1 = Flatten[Table[\[Alpha][j, n][0] == 1, {j, nMUL}, {n, nEL}]];
initCond2 = Flatten[Table[\[Lambda][i, k, q][0] == 0, {i, nMUL}, {k, nBAT}, {q, nMOD}]];

Defining a list of variables $\alpha(t), \lambda(t), \alpha^*(t), \lambda^*(t)$ to be solved for:

vars1 = Flatten[Table[\[Alpha][j, n], {j, nMUL}, {n, nEL}]];
vars2 = Flatten[Table[\[Lambda][i, k, q], {i, nMUL}, {k, nBAT}, {q, nMOD}]];
allVars = Join[vars1, vars2, vars1 /. {\[Alpha] -> \[Alpha]c}, vars2 /. {\[Lambda] -> \[Lambda]c}];

Trying to solve:

allEqsWithCond = Join[allEqs, initCond1, initCond2, initCond1 /. {\[Alpha] -> \[Alpha]c},initCond2 /. {\[Lambda] -> \[Lambda]c}];
sol1 = NDSolveValue[allEqsWithCond, allVars, {t, 0, 1}];

and getting an error:

enter image description here

Trying to solve using the Residual method:

sol2 = NDSolveValue[allEqsWithCond, allVars, {t, 0, 1}, Method -> {"EquationSimplification" -> "Residual"}];

but now getting two errors:

enter image description here enter image description here

We managed to deal with the first error by splitting variables $\alpha(t), \lambda(t), \alpha^*(t), \lambda^*(t)$ into real and imaginary parts of $\alpha(t), \lambda(t)$, however, the second error persists. We would be interested in the best approach for solving these equations.

$\endgroup$
3
  • $\begingroup$ Hm. Maybe a "homotopy method" might help. Can you connect your system by a continuous path to a simpler (maybe linear?) system that is easier to solve? Then you try to solve a sequence of systems along that path, using the solution of the previous one as initial condition until you finally arrive at your desired system (and its solution). $\endgroup$ Mar 31 '20 at 18:30
  • $\begingroup$ What is meant by the phrase connect your system by a continuous path? Difficulty of solving the presented system of equations highly depend on parameter nMUL. For nMUL = 1, the system can be reduced to non-implicit and solved, however, solution is simplistic and the true dynamics require nMUL > 1. Could this simplistic solution help us? $\endgroup$
    – JMantas
    Apr 1 '20 at 13:25
  • $\begingroup$ Yes, it could help! See, solve the equations for nMUL =1. Then increase nMUL a little and use Newton's method with the result from nMUL =1 as initial guess. When this converges, increase nMUL again a little and apply Newton's method again. Rinse and repeat. It is just not that easy to do that with NDSolve alone. You have to discretize you system; then you can use FindRoot to apply Newton's method. $\endgroup$ Apr 1 '20 at 13:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.