1
$\begingroup$

I have tried both, NDSolve and ParametricNDSolve, to tackle the following problem without success. I have looked at 2 other SE posts (here and here) that seem similar to mine, but wasn't able to resolve my problem using those. Could someone point out what I'm missing? I'd also appreciate any pointers about the deeper Wolfram Language concepts causing this issue.

The problem: I have a function f of variable x with c1 and c2 as parameters:

f[c1_,c2_,x_]:=c1^2 (1 - x c2) HeavisideTheta[c2 - x]

This function feeds the definition of the parametric model, involving an NDSolve:

model[c1_, c2_, k_] := NDSolve[
{g'[x] + (f[c1, c2, x]/k) Sin[k x + g[x]]^2 == 0, g[0] == 0}, 
g,
{x, 15/c2}]

The above NDSolve returns an InterpolatingFunction for explicit values of the arguments c1, c2 and k.

Now, the object I'm ultimately interested in is the function of k obtained by taking the last value of the InterpolatingFunction, for each value of k.

I have numeric data (Reals) in the form {{x1,y1},{x2,y2},....,{xn,yn}}. What I'd like to do is to FindFit for parameters {c1,c2} in the following sense:

FindFit[data, Last[g["ValuesOnGrid"] /. First@model[c1, c2, k]], {c1, c2}, k]

This, however, gives the error message "Endpoint 15.708/c2 in {r,15.708/c2} is not a real number". I have tried setting this problem up using ParametricNDSolve as well, but to no avail. I've attached a screen-shot of what I see.enter image description here

Improve this question
$\endgroup$
8
  • $\begingroup$ You have a syntax error in f with a bracket ) and you're using r in model - should this be x? I fixed the obvious syntax error in f but got NDSolve::underdet: There are more dependent variables, {g[r]}, than equations, so the system is underdetermined. because of the problems in your model. Also, please can you post some of your data. $\endgroup$ – flinty Jul 29 '20 at 19:35
  • $\begingroup$ You should also consider changing model to model[c1_?NumericQ, c2_?NumericQ, k_?NumericQ]:= ... $\endgroup$ – flinty Jul 29 '20 at 19:47
  • $\begingroup$ Thanks @flinty. Yes, my post had typos. I've changed the r to x and fixed the syntax error in f. I tried to make the code more readable for viewers here, and mistyped while doing that. I've also added a comment on the structure of data. $\endgroup$ – ConservedCharge Jul 29 '20 at 19:55
  • $\begingroup$ Can you provide the data, or at least a random sample? Why are you using "ValuesOnGrid" as the model? $\endgroup$ – flinty Jul 29 '20 at 19:59
  • $\begingroup$ My ultimate goal is to have a function h[k] which has c1 and c2 as parameters. For any value of k (say, k0), the way to obtain h[k0] is this: 1. Solve the ODE for g[x]. Note that this ODE for g[x] involves not only c1 and c2 as parameters, but also k0. 2. The value of h and k0, h[k0] is obtained by taking the value of the InterpolatingFunction returned by NDSolve at the rightmost point-"asymptotic" value This is the procedure for the particular value, k = k0. I want to fit this function h[k] to the data so that I can obtain what the parameters c1 and c2 in the function f[c1, c2, x] must be. $\endgroup$ – ConservedCharge Jul 29 '20 at 20:26
1
$\begingroup$

This works for me:

(* note I needed to add a Piecewise, because HeavsideTheta is non-numeric at zero *)
f[c1_, c2_, x_] := c1^2 (1 - x c2) Piecewise[{{1, c2 - x > 0}}, 0]

sol = g /. ParametricNDSolve[{g'[x] + (f[c1, c2, x]/k) Sin[k x + g[x]]^2 == 0, 
  g[0] == 0}, g, {x, 15/c2}, {c1, c2, k}];

SeedRandom[1];
data = Table[{k, 2 k^2 - RandomReal[{-2, 2}]}, {k, 0.001, 3, .1}];

(* get the endpoint value *)
getsolk[c1_?NumericQ, c2_?NumericQ, k_?NumericQ] := Last[sol[c1, c2, k]["ValuesOnGrid"]]

fit = FindFit[data, {getsolk[c1, c2, k]}, {c1, c2}, k]

(* result: {c1 -> -123.735, c2 -> -72.2024} *)
Improve this answer
$\endgroup$
1
  • $\begingroup$ Excellent! I can reproduce this at my end, and it works for my full problem as well. While I'll continue to try to understand why my approach wasn't working while this does, I'd appreciate if you had any additional insights into why my version wasn't working. $\endgroup$ – ConservedCharge Jul 29 '20 at 21:04
1
$\begingroup$

You can ask NDSolve to return the requested quantity straight away like this:

model[c1_?NumericQ, c2_?NumericQ, k_?NumericQ] := NDSolveValue[
  {g'[x] + (f[c1, c2, x]/k) Sin[k x + g[x]]^2 == 0, g[0] == 0},
  g[15/c2],
  {x, 15/c2}
];
model[1, 2, 3]

0.313396

As you can see, this will return the value of g at the point 15/c2 instead of the full interpolation function.

Improve this answer
$\endgroup$
1
  • $\begingroup$ Thanks for the tip, @Sjoerd; this definitely makes it cleaner. $\endgroup$ – ConservedCharge Jul 30 '20 at 5:13

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.