1
$\begingroup$

I am a self-taught beginner trying to use Mathematica for the first time. If you wouldn't mind, I would like to ask for help with the code I am working on as I keep running into multiple issues when using NDSolve to solve for a problem where given four ODEs and some terminal conditions, I am asked for the time (tf) at which the terminal conditions should be reached.

The equations are somewhat hypothetical for tractability of our real equations that deal with water, so there's not much of a meaningful story behind them.

The issues I face are: Depending on the choice on the FindRoot option for the starting point (t0 + x0), I end up with NDSolve::ndsz, NDSolve::ndsv:, NDSolve::mxst, etc.

Also, I believe that the equations are stiff, and should have a singularity at W=0.

The main problem, however, is that despite all attempts, I still couldn't get any value for tf to solve the problem

With these, I was wondering if anyone could help. Your input would be greatly appreciated!

(*Parameters and Functions*)
t0 = 1991; r = .03; pb = 10; W0 = 5000; s = .05; w = 5000; cq = 3; alpha = 1; g = .02; eta = .3; cd = .947

l = 0.00005*W[t]^2
lprimew = D[l, W[t]]
evap = -.00004 W[t]^2
evapprimew = D[evap, W[t]]
ct = 10*S[t]
ctprime = D[ct, S[t]]
damage = 10*S[t]^2
damageprime = D[damage, S[t]]
q = alpha*E^(g (t - t0)) (p[t] + cd)^(-eta)

(*Terminal W, S, Mu*)
Wprox = W[t] /. Solve[w - l - evap == 0, W[t]]
Wmax = First[Select[W[t] /. Solve[w - l - evap == 0, W[t]], # > 0 &]]

Optans = Select[{\[Mu][t], W[t], S[t]} /. Solve[{pb - cq - ct == -1/(evapprimew + lprimew + r)*((S[t]*damageprime/W[t]) + ((s*w - S[t]*l)/W[t])*ctprime + (\[Mu][t]/(W[t])^2)*(2*s*w - S[t]*(w - evap + evapprimew*W[t]))), \[Mu][t] == -(ctprime*(w - evap - l) + damageprime)/(r + (w - evap)/W[t]), S[t] == s*w/(w - evap - 2*l)}, {\[Mu][t], W[t], S[t]}], Re[#[[1]]] < 0 && Wmax >=  Re[#[[2]]] > 0 && 1 >= Re[#[[3]]] > 0 &]
Wopt = Optans[[1, 2]]
Sopt = Optans[[1, 3]]
\[Mu]opt = Optans[[1, 1]]

(*Differential Equations*)
peqn = p'[t] == (p[t] - cq - ct)*(evapprimew + lprimew + r) + (((S[t]*damageprime/W[t]) + ((s*w - S[t]*l)/W[t])*ctprime + (\[Mu][t]/(W[t])^2)*(2*s*w - S[t]*(w - evap + evapprimew*W[t]))));
Weqn = W'[t] == w - q - evap - l;
Seqn = S'[t] == (s*w - S[t]*(q + l))/W[t];
mueqn = \[Mu]'[t] == r*\[Mu][t] + ctprime*q + damageprime + \[Mu][t]*(q + l)/W[t];

(*Numerical Solution*)
system[x_] := {peqn, Weqn, Seqn, mueqn, p[x] == pb, W[x] == Wopt, S[x] == Sopt, \[Mu][x] == \[Mu]opt}
solution[tf_] := NDSolve[system[tf], {p[t], W[t], S[t], \[Mu][t]}, {t, t0, tf}, Method -> {"StiffnessSwitching", "NonstiffTest" -> False}, AccuracyGoal -> 4, PrecisionGoal -> 4]
Winitial[tf_?NumberQ] := First[(W[t] /. solution[tf]) /. t -> t0]
tfopt = tf /. FindRoot[Winitial[tf] == W0, {tf, t0 + 8}]

Result = solution[tfopt];
Plot[Evaluate[W[t] /. Result], {t, t0, tfopt}]
$\endgroup$
  • $\begingroup$ Welcome to Mathematica.SE! I hope you will become a regular contributor. To get started, 1) take the introductory tour now, 2) when you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge, 3) remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign, and 4) give help too, by answering questions in your areas of expertise. $\endgroup$ – bbgodfrey Dec 30 '15 at 18:20
1
$\begingroup$

A Plot of Winitial

Plot[Winitial[tf], {tf, t0, t0 + 2}]

enter image description here

shows that it is real over only a narrow range of tf, and FindRoot in combination with NDSolve has difficulty outside that range. So, try

tfopt = tf /. FindRoot[Winitial[tf] == W0, {tf, 1991.5}]
(* 1991.73 *)

which is, I believe, the correct answer.

$\endgroup$
  • $\begingroup$ Thank you so much!! I tried initially higher values as starting points (such as t0 + 9), based on my computations on optimal q, and then moved in multiples of 5 (tried 4). I didn't expect it to be such a small number! I just want to ask how you thought of such a small starting number to try for the FindRoot command. Anyway, thank you again!!! $\endgroup$ – Anonymous Dec 31 '15 at 1:54
  • $\begingroup$ @Anonymous I read a guess off the plot. $\endgroup$ – bbgodfrey Dec 31 '15 at 4:00

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.