2
$\begingroup$

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?

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

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"]}

Figure 1

$\endgroup$

Your Answer

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

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