I would like to fit a set of data (measured value vs time) with a model that contains a differential equation and an integral that has to be evaluated numericaly. The integral also contains Bessel functions (first and second kind) evaluated close to zero which I think does not help to solve the problem because of divergence.

The model is the following: $$ \frac{dT\!\left(t\right)}{dt}=-\frac{16 c_0 Q}{\pi}\int_0^\infty \frac{\exp\left(-Qu^2t\right)}{u \left(J_0^2\left(u\rho\right)+Y_0^2\left(u\rho\right)\right)}T\!\left(t\right)^2 du $$ $u$ is just a dummy integration variable, and the parameters to be optimised are $c_0, \rho$ and $Q$.

Here is the piece of code that I wrote, however without success:

(* Fit model *) 
model = ParametricNDSolveValue[{T'[t] ==
                                NIntegrate[-((16 c Q)/Pi) Exp[-Q u^2 t]/
                                           (BesselJ[0,u p]^2 + BesselY[0, u p]^2) T[t]^2,
                               {u, 0, ∞}], T[0] = 1}, T, {t, 0, 2*^-6}, {c, Q, p}]
(* fitting *) 
fit = NonlinearModelFit[data,model[c, Q, p][t], {{c, 1*^15}, {Q, 1*^-6}, {p, 2*^-9}}, t]

There are many errors, the first being that NIntegrate is evaluated at a non-numerical value.

Do you have suggestions to improve the code?

Here is a test data sample:

data = {{0, 1.05634}, {1/10000000, 0.752195}, {1/5000000, 0.495366},
        {3/10000000, 0.484146}, {1/2500000, 0.342561}, {1/2000000, 0.43622},
        {3/5000000, 0.235122}, {7/10000000, 0.128049}, {7.95*10^-7, 0.147317},
        {9/10000000, 0.192805}, {1/1000000, 0.14561}, {1.1*10^-6, 0.32939},
        {1.2*10^-6, 0.228659}, {1.3*10^-6, -0.0340244}, {1.4*10^-6, 0.0121951},
        {1.5*10^-6, 0.131707}, {1.6*10^-6, 0.074878}, {1.7*10^-6, 0.30939},
        {1.8*10^-6, 0.122561}, {1.9*10^-6, -0.00195122}, {2.005*10^-6, 0.0857317}};

**** EDIT ****

I have continued working on this problem. Now I first solve the integral, and I pass it to ParametricNDSolve. I do not think that evaluating Bessel functions close to 0 is a problem as the evaluation of the integral alone works well. The code is thus as follow:

f[c_?NumericQ, Q_?NumericQ, p_?NumericQ, t_?NumericQ] := 
  NIntegrate[- ((16 c Q)/Pi) Exp[-Q u^2 t]/(BesselJ[0, u p]^2 + BesselY[0, u p]^2), 
 {u, 0, ∞}];

model = T /.ParametricNDSolve[{T'[t] == f[c, Q, p][t] T[t]^2, T[0] == 1}, 
T, {t, 0, 100}, {c, Q, p}];    

fit = FindFit[data[[All, 2]],model[c, Q, p][t], {{c, 1*^15}, {Q, 1*^-6}, {p, 2*^-9}}, t];

Now I have the following error message: "ParametricNDSolve::pdvar Dependent variables {T,f[c,Q,p]} cannot depend on parameters {c,Q,p}". I find it surprising because in principle one wants the dependent variables to depend on some parameters.

  • $\begingroup$ It is because the result of function T inside the integral is unknown at the moment of integral evaluation. Try to write your own Runge-Kutta DE solving algorithm as a function of c, q and P. The first step of algorithm is assigning initial value for T function. Thus the result of T[0] will not be non-numerical value. $\endgroup$ – Sergey Oct 3 '17 at 5:39
  • $\begingroup$ Integro-differential equations are not yet supported, I believe. $\endgroup$ – J. M. is away Oct 3 '17 at 5:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.