# [Solved] Problem using Piecewise and Module inside NDSolve.

Posted 9 months ago
1067 Views
|
5 Replies
|
4 Total Likes
|
 (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?
5 Replies
Sort By:
Posted 9 months ago
 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 9 months ago
 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 9 months ago
 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