0

I'm trying to find the root (po) of a complicated problem (where s is a parameter):

First, I need to define some functions, which involve numerically solving a differential equation.

Needs["DifferentialEquations`NDSolveProblems`"];
Needs["DifferentialEquations`NDSolveUtilities`"];
Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"];

a[po_, s_] := NDSolve[{3*F[p] + p*F'[p] - 2*s*(3*F''[p] + p*F'''[p]) == 0, F[po] == 1, F'[po] == 0, 
F''[po] == (-1 + 3*po - 3*po^2 - 2*s)/(2*s*po*(1 - 2*po)), WhenEvent[F[p] == 0, "StopIntegration"]}, F, {p, 0.01}]

plownum[po_, s_] := InterpolatingFunctionDomain[First[F /. a[po, s]]][[1, 1]]

plowdensity[po_, s_] := First[F'[plownum[po, s]] /. a[po, s]]

plowprofit[po_, s_] := plownum[po, s]^2 (1 + 2 s*plowdensity[po, s])

rootpo[po_, s_] := plowprofit[po, s] - (po^2*(1 - 2*po + po^2 + 2*s))/(2 - 4*po)

Without explaining in detail what I'm trying to achieve, everything works perfectly up to that point (since I'm not extremely familiar with the intricacies of Mathematica, even that took me a while). The function rootpo[po,s] evaluates very quickly and without problems.

However, when trying to find the value of po where rootpo[po,s] == 0, I get a problem when using FindRoot:

poeq[s_] := FindRoot[rootpo[po, s], {po, 0.1}]

poeq[0.2]

leads to a large error output:

NDSolve::ndsv: Cannot find starting value for the variable F.
ReplaceAll::reps: {NDSolve[{...},F,{p,0.01}]} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing.

etc.

I tried different things (including looking up solutions to similar problems and adjusting the code), but unfortunately I still couldn't get the FindRoot[.] to run.

Can anybody help me? I'm also happy about suggestions how to clean up the code.

Many thanks in advance!

Martin
  • 101
  • 3
  • 1
    You might want some ?NumericQ on the function arguments in the definitions, like a[po_?NumericQ, s_?NumericQ] := .... – Michael E2 Jun 22 '17 at 19:16
  • Thanks! I just tried that (changed the arguments in every function definition), but unfortunately I get the same error message. – Martin Jun 22 '17 at 19:23
  • you should update the question with that fix. Be sure to restart your kernel after making such change as well. – george2079 Jun 22 '17 at 19:26
  • Also your NDSolve domain doesn't look right, perhaps should be {p,po,upperlim} where the upper limit is some sufficiently large number or Infinity presuming you want it to just go until WhenEvent stops it. ? – george2079 Jun 22 '17 at 19:37
  • taking back the last comment, that seems to work, NDSolve evidently takes a limit from your initial condition, so you implicitly have {p, 0.01, po}. I cant see where that feature is documented though. – george2079 Jun 22 '17 at 19:54
  • I didn't restart the kernel initially, and it seems that things do change quite a bit after the suggested modification. I will look into that tomorrow. Regarding your other comments: The NDSolve seems to work just fine, even though I might be lucky there. Thanks again! – Martin Jun 22 '17 at 22:40
  • 2
    Martin, I should've inserted a ClearAll[a] to wipe out the previous definition: ClearAll[a]; a[po_?NumericQ,.... Restarting the kernel clears all definitions, too, of course. -- @george2079 {x, x1} is an old legacy form that implies {x, x0, x1} where x0 is taken from the IC (so {p, po, 0.01} implied in the OP's code). – Michael E2 Jun 22 '17 at 23:25
  • Thanks a bunch - adding ?NumericQ to the function definitions did work indeed (after restarting the kernel). I thus consider the problem as solved... – Martin Jun 26 '17 at 11:14

0 Answers0