Why can't I solve a system of four third-order partial differential equations using NDSolve?

I solve the problem of optimal control. First, I solve a system of two equations (with respect to the variables p, q, V, n) with already optimal conditions (Kopt and Sopt) to check whether it is possible to solve such a system at all. It turns out that everything is OK, but the non-linear part of the equation cannot be discarded, it is important. I am using the NDSolve function with the StiffnessSwitching specification:

Kopt[t_] :=
6.03626835*10^17*t^7 - 4.51848425*10^16 *t^6 +
1.36700872*10^15*t^5 - 2.13240902*10^13*t^4 +
1.80975334*10^11*t^3 - 7.99509911*10^8*t^2 +
1.51825005*10^6 *t^1 + 1.87782986*t^0;

Sopt[t_] :=
0. - 90989.2 t + 1.39914*10^9 t^2 - 8.76585*10^12 t^3 +
2.99045*10^16 t^4 - 6.32825*10^19 t^5 + 9.0249*10^22 t^6 -
9.14249*10^25 t^7 + 6.81239*10^28 t^8 - 3.82682*10^31 t^9 +
1.65063*10^34 t^10 - 5.54261*10^36 t^11 + 1.46314*10^39 t^12 -
3.0549*10^41 t^13 + 5.05619*10^43 t^14 - 6.62179*10^45 t^15 +
6.81692*10^47 t^16 - 5.44703*10^49 t^17 + 3.30728*10^51 t^18 -
1.47379*10^53 t^19 + 4.54318*10^54 t^20 - 8.65281*10^55 t^21 +
7.66886*10^56 t^22;

eq1 = D[n[x, t], t] + D[n[x, t]*V[x, t], x] == 0;
eq2 = D[V[x, t], t] + Kopt[t]*x -
2*(2*250^2)*Sopt[t]*Sin[k*x]*Cos[k*x]*k + V[x, t]*D[V[x, t], x] +
g*D[n[x, t], x] -
D[ D[Sqrt[n[x, t]], x, x]/(2*Sqrt[n[x, t]]), x] == 0;

eqo1 = D[n[x, t], t] + D[n[x, t]*V[x, t], x] == 0;
eqo2nl = D[V[x, t], t] + Kopt[t]*x -
8*(2*250^2)*Sopt[t]*Sin[k*x]*Cos[k*x]*k + V[x, t]*D[V[x, t], x] +
g*D[n[x, t],
x] - ((D[n[x, t], x])^3 -
2*n[x, t]*D[n[x, t], x]*D[n[x, t], {x, 2}] +
D[n[x, t], {x, 3}]*(n[x, t])^2)/(4*(n[x, t])^3) == 0;

cond1 = n[x, t] == ((Sqrt[Kzero]/Pi)^(1/2))*
Exp[-0.5*Sqrt[Kzero]*x^2] /. t -> 0;
cond2 =  n[x, t] == 0. /. x -> -10;
cond3 =  n[x, t] == 0. /. x -> 10;

cond4 = V[x, t] == 0. /. t -> 0;
cond5 = V[x, t] == 0. /. x -> -10;
cond6 = V[x, t] == 0. /. x -> 10;

sol5 = NDSolve[{eqo1, eqo2nl, cond1, cond2, cond3, cond4, cond5,
cond6}, {V, n}, {x, -10, 10}, {t, 0, Tv},
Method -> {"StiffnessSwitching",
Method -> {"ExplicitRungeKutta", Automatic}}, AccuracyGoal -> 1,
PrecisionGoal -> 1];


Next, I solve four equations in order to find the control parameter K (depends on q), while Sopt is still known. The initial conditions are the same and as banal as possible. Add two more equations and K[t] formula:

K[t_] := -1*Integrate[q[x, t]*x^2, {x, -10, 10}];

eqo3 = D[q[x, t], t] == - n[x, t]*D[p[x, t], x]  -
V[x, t]*D[q[x, t], x];

eqo4 = D[p[x, t], t] +  D[p[x, t], x]*V[x, t] + g*D[q[x, t], x] == 0;


The linearized system has a solution, but not a nonlinear one! There are many problems, including the order of the equation:

The spatial derivative order of the PDE may not exceed two.


I try to remove the third order derivative from the nonlinear part:

The maximum derivative order of the nonlinear PDE coefficients for \
the Finite Element Method is larger than 1. It may help to rewrite \
the PDE in inactive form.


Why does he even use only FEM now? Why was a nonlinear system of two equations solved normally?

• Parameters k, g, Tv, Kzero, Sopt, Kopt not defined. Feb 10, 2020 at 22:49
• k = (20.3/589)*10^3, g = 1039, Kzero = 10. Kopt and Sopt are interpolation polynomials. But they do not bring the main problem. Feb 11, 2020 at 7:50
• Add to your post this data including functions Kopt, Sopt. Feb 11, 2020 at 12:45
• Did it, thanks. Feb 11, 2020 at 14:17
• Ok! What is there Tv? Feb 11, 2020 at 17:26

We use Knopt[t] instead of K[t] (K is a system symbol) and initial data p[x, 0] == Cos[Pi x/20], q[x, 0] == 0 instead of c7 = p[x, t] == 0 /. t -> Tv; c10 = q[x, t] == 0 /. t -> Tv. Then the numerical solution converges:

Sopt[t_] :=
0. - 90989.2 t + 1.39914*10^9 t^2 - 8.76585*10^12 t^3 +
2.99045*10^16 t^4 - 6.32825*10^19 t^5 + 9.0249*10^22 t^6 -
9.14249*10^25 t^7 + 6.81239*10^28 t^8 - 3.82682*10^31 t^9 +
1.65063*10^34 t^10 - 5.54261*10^36 t^11 + 1.46314*10^39 t^12 -
3.0549*10^41 t^13 + 5.05619*10^43 t^14 - 6.62179*10^45 t^15 +
6.81692*10^47 t^16 - 5.44703*10^49 t^17 + 3.30728*10^51 t^18 -
1.47379*10^53 t^19 + 4.54318*10^54 t^20 - 8.65281*10^55 t^21 +
7.66886*10^56 t^22;
Knopt[t_] := -Integrate[q[x, t]*x^2, {x, -10, 10}];
k = (20.3/589)*10^3; g = 1039; Kzero = 10.; Tv = 3.9*10^(-6);
eq1 = D[n[x, t], t] + D[n[x, t]*V[x, t], x] == 0;
eq2 = D[V[x, t], t] + Kopt[t]*x -
2*(2*250^2)*Sopt[t] Sin[k*x]*Cos[k*x]*k + V[x, t]*D[V[x, t], x] +
g*D[n[x, t], x] -
D[D[Sqrt[n[x, t]], x, x]/(2*Sqrt[n[x, t]]), x] == 0;

eqo1 = D[n[x, t], t] + D[n[x, t]*V[x, t], x] == 0;
eqo2nl = D[V[x, t], t] + Kopt[t]*x -
8*(2*250^2)*Sopt[t]*Sin[k*x]*Cos[k*x]*k + V[x, t]*D[V[x, t], x] +
g*D[n[x, t],
x] - ((D[n[x, t], x])^3 -
2*n[x, t]*D[n[x, t], x]*D[n[x, t], {x, 2}] +
D[n[x, t], {x, 3}]*(n[x, t])^2)/(4*(n[x, t])^3) == 0;

cond1 = n[x, t] == ((Sqrt[Kzero]/Pi)^(1/2))*
Exp[-0.5*Sqrt[Kzero]*x^2] /. t -> 0;
cond2 = n[x, t] == 0. /. x -> -10;
cond3 = n[x, t] == 0. /. x -> 10;

cond4 = V[x, t] == 0. /. t -> 0;
cond5 = V[x, t] == 0. /. x -> -10;
cond6 = V[x, t] == 0. /. x -> 10;

eqo3 = D[q[x, t], t] == -n[x, t]*D[p[x, t], x] - V[x, t]*D[q[x, t], x];

eqo4 = D[p[x, t], t] + D[p[x, t], x]*V[x, t] + g*D[q[x, t], x] ==
0; eqo2n2 =
D[V[x, t], t] + Knopt[t]*x -
8*(2*250^2)*Sopt[t]*Sin[k*x]*Cos[k*x]*k + V[x, t]*D[V[x, t], x] +
g*D[n[x, t],
x] - ((D[n[x, t], x])^3 -
2*n[x, t]*D[n[x, t], x]*D[n[x, t], {x, 2}] +
D[n[x, t], {x, 3}]*(n[x, t])^2)/(4*(n[x, t])^3) == 0;

ic = {n[x, 0] == ((Sqrt[Kzero]/Pi)^(1/2))*Exp[-0.5*Sqrt[Kzero]*x^2],
V[x, 0] == 0, p[x, 0] == Cos[Pi x/20], q[x, 0] == 0};
bc = {cond2, cond3, cond5, cond6, p[-10, t] == 0, p[10, t] == 0,
q[-10, t] == 0, q[10, t] == 0};

sol5 = NDSolve[{eqo1, eqo2n2, eqo3, eqo4, ic, bc}, {V, n, p,
q}, {x, -10, 10}, {t, 0, Tv}];
{DensityPlot[V[x, t] /. sol5, {x, -10, 10}, {t, 0, Tv},
ColorFunction -> "Rainbow", PlotLegends -> Automatic,
PlotLabel -> "V"],
DensityPlot[n[x, t] /. sol5, {x, -10, 10}, {t, 0, Tv},
ColorFunction -> "Rainbow", PlotLegends -> Automatic,
PlotRange -> All, PlotLabel -> "n"],
DensityPlot[p[x, t] /. sol5, {x, -10, 10}, {t, 0, Tv},
ColorFunction -> "Rainbow", PlotLegends -> Automatic,
PlotLabel -> "p"],
DensityPlot[q[x, t] /. sol5, {x, -10, 10}, {t, 0, Tv},
ColorFunction -> "Rainbow", PlotLegends -> Automatic,
PlotRange -> All, PlotLabel -> "q"]}