Message Boards Message Boards

0
|
3082 Views
|
5 Replies
|
4 Total Likes
View groups...
Share
Share this post:

[Solved] Problem using Piecewise and Module inside NDSolve.

Marco Maria Molinari
Posted 4 years ago

(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?
POSTED BY: Marco Maria Molinari
5 Replies
Sort By:

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
POSTED BY: Neil Singer

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
POSTED BY: Marco Maria Molinari

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
POSTED BY: Neil Singer

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
POSTED BY: Marco Maria Molinari

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
POSTED BY: Neil Singer
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Tag limit exceeded
Note: Only the first five people you tag will receive an email notification; the other tagged names will appear as links to their profiles.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract