1
$\begingroup$

Observed an error when calling a Module from FindRoot.

I created a module which first solves a differential equation and then returns the value of that differential equation when the independent variable, time, has a value of 1800 (seconds). The model accepts an input which is one of the parameters of the system, vSys.

Note that this example is just the first step of creating a more complete system. The goal is to provide some tools to find the optimal parameters of the system given various constraints. So, in a later implementation would like to pass additional parameters, such as other system parameters and the target time. Also, in later implementations the differential equation will not have an analytic solution.

Below is the module:

p30[vSysInput_] := 
 Module[{eq2A, p, t, rGas, Tgas, vSys, parms, sol},
    eq2A = p'[t] == rGas Tgas  / vSys (-0.64 *0.0012 p[t] );
    parms = {rGas -> 8.3144, Tgas -> 273, vSys -> vSysInput};
    sol = NDSolve[{eq2A /. parms, p[0] == 100000}, p[t], {t, 0, 1800}];
    p[t] /. sol[[1]] /. t -> 30 * 60
  ]

In some contexts the module appears to work as expected. For example, the module can be called to determine a value for a given input and the module can be called from Plot.

{#, p30[#]} & /@ {800, 1000, 1200} // TableForm]

Mathematica graphics

Plot[p30[vol], {vol, 500, 1500}]

Mathematica graphics

But problems appear when the module is used with FindRoot:

FindRoot[ p30[vol] - 4000, {vol, 800, 1200}]

Mathematica graphics

The problem appear to be related to the parameter passed to the model remaining a symbol rather than changing to a numeric value.

Using _NumericQ, resulted in a different set of problems.
Mathematica graphics

Mathematica graphics

Thanks for your help. Thanks everyone -- works much better with _?, rather than just _

$\endgroup$
  • $\begingroup$ Your images do not show correctly. If you are using the SETools palette, use the "Image" button, not the "Selected Notebook" one. $\endgroup$ – MarcoB May 18 '15 at 19:10
  • 1
    $\begingroup$ Try Clear[p30]; p30[vSysInput_?NumericQ] := ... . $\endgroup$ – ilian May 18 '15 at 19:29
3
$\begingroup$

I think in the latest revision there may be missing a PatternTest (or ?):

In[4]:= Clear[p30]; 
        p30[vSysInput_?NumberQ] := 
          Module[{eq2A, p, t, rGas, Tgas, vSys, parms, sol}, 
            eq2A = p'[t] == rGas Tgas/vSys (-0.64*0.0012 p[t]);
            parms = {rGas -> 8.3144, Tgas -> 273, vSys -> vSysInput};
            sol = NDSolve[{eq2A /. parms, p[0] == 100000}, p[t], {t, 0, 1800}];
            p[t] /. sol[[1]] /. t -> 30*60]
        FindRoot[p30[vol] - 4000, {vol, 800, 1200}]

Out[5]= {vol -> 974.817}

It's a somewhat frequently asked question, please see this writeup.

$\endgroup$

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.