[Solved] Problem using Piecewise and Module inside NDSolve.

Questions
, ,
1

(Index i[z])
j = i[z[h]];

 (*ISA model*)
 out = T0[[j]] + Lstd[[j]]/1000 (z[h] - zstd[[j]]);
 
 Return[out];
 ]
Pstd[h_] := Module[
  {Lstd, zstd, T0, pb, j, out},
  
  (*Coefficients definitions*)
  Lstd = {-6.5, 0., 1., 2.8, 0, -2.8, 0};
  zstd = {0., 11000., 20000., 32000., 47000., 51000, 80000};
  T0 = {288.15, 216.65, 216.65, 228.65, 270.65, 270.65, 192.247};
  pb = {101325, 22694, 5544.63, 900.192, 120.096, 76.3148, 1.17443};
  
  (*Index i[z]*)
  j = i[z[h]];
  
  (*ISA model*)
  out = If[Lstd[[j]] == 0, 
    pb[[j]] Exp[-g0/(R*T0[[j]]) (z[h] - zstd[[j]])], 
    pb[[j]] (1 + Lstd[[j]]/(T0[[j]]*1000) (z[h] - zstd[[j]]))^(-(g0/(
      Lstd[[j]]*0.287)))];
  Return[out];
  ]
\[Rho]std[h_] := Module[
  {out},
  out = Pstd[h]/(R*Tstd[h]);
  
  Return[out];
  ]
NDSolve[{V'[t] == -0.5* \[Rho]std[h[t]]*(V[t])^2 ((Cd*S)/m) + 
    g[h[t]] Cos[\[Theta][t]],
  \[Theta]'[t] == -0.5*\[Rho]std[h[t]]*V[t] ((Cl*S)/m) + (
    g[h[t]]*Sin[\[Theta][t]])/V[t] - V[t]/(Rt + h[t]) Cos[\[Theta][t]],
  h'[t] == -V[t]*Sin[\[Theta][t]],
  V[0] == 11140,
  \[Theta][0] == 6.92,
  h[0] == 120*10^3
  }, {V, \[Theta], h}, {t, 0, 7200}]

I know the problem is that I can’t provide a numerical value to \[Rho]std[h[t]] inside Ndsolve. Advices?

2

Marco,

Your code is not complete so I can’t say why you are getting an error, however, I took a sample from Mathematica’s documentation and added a piecewise function in a similar way and this works fine. I used Round[t] to create an index into a list because the code that finds your index is missing. func[] is a discontinuous function that looks up values in a table. newX[t] is a module with an If condition and I use them in NDSolve.

func[t_] := 
 Module[{tab = {1.1, 2.2, 3.3, 4.4, 5.5}, temp}, temp = Round[t]; 
  tab[[temp]]]

newX[t_] := 
 Module[{out, tab }, out = If[x[t] >= 0, x[t], func[t]*x[t]]; 
  Return[out]]

NDSolve[{Derivative[1][x][t] == -3 (newX[t] - y[t]), 
    Derivative[1][y][t] == -x[t] z[t] + 26.5 x[t] - y[t], 
    Derivative[1][z][t] == x[t] y[t] - z[t], x[0] == z[0] == 0, 
    y[0] == 1}, {x, y, z}, {t, 0, 5}];

Plot[Evaluate[{x[t], y[t], z[t]} /. First[%]], {t, 0, 5}]

Your bug must be something with the way you are doing your index. To be more helpful, I would need to see the rest of the code that is cut off at the top.

Regards,

Neil

3

Dear Neil,

First of all thanks for the reply and sorry if I cut off the top of the script, my bad. I have now entered the whole script. I probably think that the bug is due to the fact that the rho variable as it is defined does not depend directly on the time. Even seeing your script that has been helpful to me, I am not able to fix it.

g0 = 196133/20000;
R = 287.05;
\[Gamma] = 1.4;
Cpmax = ((\[Gamma] + 1)^2/(4 \[Gamma]))^(\[Gamma]/(\[Gamma] - 
    1))*4/(\[Gamma] + 1);
Rt = 6325766; 

g[z_] := g0*(Rt/(Rt + z))^2;
Cl = 0;
Cd = 1.67929;
S = 12;
m = 5424.9;

z[h_] := r0*h/(r0 + h) /. r0 -> 6325766.
i[z_] := Piecewise[{{1, z <= 11000}, {2, 11000 < z <= 20000}, {3, 
     20000 < z <= 32000}, {4, 32000 < z <= 47000}, {5, 
     47000 < z <= 51000}, {6, 51000 < z <= 80000}, {7, 
     80000 < z < 120000}}];
Tstd[h_] := Module[
  {Lstd, zstd, T0, j, out},
  
  (*Coefficients definitions*)
  Lstd = {-6.5, 0., 1., 2.8, 0, -2.8, 0};
  zstd = {0., 11000., 20000., 32000., 47000., 51000, 80000};
  T0 = {288.15, 216.65, 216.65, 228.65, 270.65, 270.65, 192.247};
  
  (*Index i[z]*)
  j = i[z[h]];
  
  (*ISA model*)
  out = T0[[j]] + Lstd[[j]]/1000 (z[h] - zstd[[j]]);
  
  Return[out];
  ]
Pstd[h_] := Module[
  {Lstd, zstd, T0, pb, j, out},
  
  (*Coefficients definitions*)
  Lstd = {-6.5, 0., 1., 2.8, 0, -2.8, 0};
  zstd = {0., 11000., 20000., 32000., 47000., 51000, 80000};
  T0 = {288.15, 216.65, 216.65, 228.65, 270.65, 270.65, 192.247};
  pb = {101325, 22694, 5544.63, 900.192, 120.096, 76.3148, 1.17443};
  
  (*Index i[z]*)
  j = i[z[h]];
  
  (*ISA model*)
  out = If[Lstd[[j]] == 0, 
    pb[[j]] Exp[-g0/(R*T0[[j]]) (z[h] - zstd[[j]])], 
    pb[[j]] (1 + Lstd[[j]]/(T0[[j]]*1000) (z[h] - zstd[[j]]))^(-(g0/(
      Lstd[[j]]*0.287)))];
  Return[out];
  ]
\[Rho]std[h_] := Module[
  {out},
  out = Pstd[h]/(R*Tstd[h]);
  
  Return[out];
  ]

sol = NDSolve[{V'[t] == -0.5* \[Rho]std[h[t]]*(V[t])^2 ((Cd*S)/m) + 
     g[h[t]] Sin[\[Theta][t]],
   \[Theta]'[
     t] == -0.5*\[Rho]std[h[t]]*
      V[t] ((Cl*S)/m) + (g[h[t]]*Cos[\[Theta][t]])/
      V[t] - (V[t]/(Rt + h[t])) Cos[\[Theta][t]],
   h'[t] == -V[t]*Sin[\[Theta][t]],
   V[0] == 11140,
   \[Theta][0] == 6.92*\[Pi]/180,
   h[0] == 120000
   }, {V, \[Theta], h}, {t, 0, 7200}]

Kind regards,

Marco

4

Marco,

The problem is that your [Rho]std function uses Part deep inside and it will fail on any non numeric argument (because it can’t take Part of the expression until it is evaluated numerically). The easiest way to deal with this is to “protect” your function from being evaluated with non-numerical expressions:

\[Rho]std[h_Real] := Module[{out}, out = Pstd[h]/(R*Tstd[h]);
  Return[out];]

This change will get you an answer because the function will not match a pattern and evaluate until the argument to it is Real.

Note: This will fail on Integers, so other ways to protect functions: (Real only, Integer Only, Any Number):

funR[x_Real] := x
funI[x_Integer] := x
funN[x_ /; NumericQ[x]] := x

I only needed to use Real above because NDSolve and the other Numerical functions always use Reals.

Regards,

Neil

5

Mr Singer,

thank you very much for the help, now the script works but above all you have expanded my knowledge on “protection” and “conditions” that it is necessary to impose on functions. I really appreciate it.

Thanks again.

Kind regards,

Marco

6

Marco,

Glad I could help!

You will see this more often with NSolve because NSolve first tries to solve an equation analytically and then does a numerical solution.

Regards,

Neil