using FindRoot with a Module that takes a parameter

Question

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]

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

But problems appear when the module is used with FindRoot:

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

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.

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

Answer 1

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.