How to correctly write a PDE solving module for NonlinearModelFit

I need to fit a model defined by numerical operations to test data and find the parameters of the PDE. The problem is that I have to NIntegrate the result obtained by NDSolve. I prepared the following module (using the simplest PDE):

model[M0_?NumberQ, k_?NumberQ, dq_?NumberQ, tmax_?NumericQ] :=
Module[{funLoad, sol, spInit, sp, z, t, H = 0.020, gw = 10000,
solSP, solSett},
spInit[z_] := 0;
funLoad[t_] := If[t >= 0 , -dq, 0];
sol := Quiet@
NDSolve[{k/gw D[sp[z, t], {z, 2}] == (1/M0) D[sp[z, t], t],
sp[z, 0] == spInit[z], sp[0, t] == funLoad[t],
sp[H, t] == funLoad[t]}, sp, {z, 0, H}, {t, 0, tmax},
Method -> "BDF"];
solSP := First@Evaluate[sp /. sol];
solSett[tt_?NumericQ] := (-1/M0) *
NIntegrate[solSP[z, tt], {z, 0, H}];
solSett
]


Example data points:

data = {{8102., 0.00029}, {16294., 0.00034}, {24486.,
0.000379}, {32678., 0.000416}, {40870., 0.000437}, {49062.,
0.00045000000000000004}, {57254., 0.000471}, {66678.93400000001,
0.000486}, {87158.93400000001, 0.000512}};


I can test the module - it produces good results:

tmax = First@Last@data;
model[1 10^6, 10^-10, 100000, tmax][50]
(*0.000154844*)
Plot[model[2 10^6, 10^-12, 100000, tmax][t], {t, 0, tmax},
Epilog :> Point[data], PlotRange -> All]


So now I'd like to use NonlinearModelFit to find the parameters k and M0. I limited the number of steps to only 2 because the calculation is quite slow. Also, I added StepMonitor to see how the parameters change during fitting.

nlm = NonlinearModelFit[
data, {model[C1, C2, 100000, tmax][t],
C1 > 1 10^6 && C1 < 10 10^6 && C2 > 10^-14 && C2 < 10^-10}, {C1,
C2}, t, Method -> Automatic, MaxIterations -> 2,
StepMonitor :> Print["C1 = ", C1, "    C2 = ", C2]]
(*C1 = 1.*10^6    C2 = 9.0001*10^-11*)
(*C1 = 1.*10^6    C2 = 9.0001*10^-11*)


Unfortunately, I get:

nlm["BestFitParameters"]
(*{C1 -> 1.*10^6, C2 -> 1.*10^6}*)


So the second parameters is outside the range specified by constraints. Also, if I calculate more steps there is no change in C1 and C2.

Clearly, I made a mistake in writing the module model. The question is: where are my mistakes. Also I wonder: what can I do to speed up the calculations? NIntegrate takes quite a lot of time.

I think the issue concerns inadequate precision in large part because the two parameters (C1 and C2) are orders of magnitude apart. So the first step is to get the two parameters on the same scale. Also, scaling the dependent variable wouldn't hurt either.

I suspect that the function model could also use some numerical stability/precision improvements but I don't want to touch that function. Below I just make NonlinearModelFit a bit more numerically stable:

(* Scale the dependent variable *)
data2 = data;
data2[[All, 2]] = 10^4 data[[All, 2]];

(* Run regression with scaled coefficients *)
nlm = NonlinearModelFit[
data2, {10^4 model[C1 10^6, C2 10^(-11), 100000, tmax][t],
1 < C1 < 10 && 10^(-3) < C2 < 10}, {{C1, 4.43 }, {C2, 0.977}}, t,
StepMonitor :> Print["C1 = ", C1, "    C2 = ", C2],
MaxIterations -> 5];

C1 = 4.43       C2 = 0.977
C1 = 4.42592    C2 = 0.977203
C1 = 4.43656    C2 = 0.974895
C1 = 4.4366     C2 = 0.974945
C1 = 4.43649    C2 = 0.968321
C1 = 4.43594    C2 = 0.968415

(* Summary statistics *)
nlm["EstimatedVariance"]
(* 0.145749 *)
nlm["BestFitParameters"]
(* {C1 -> 4.43594, C2 -> 0.968415} *)

(* Make predictions *)
{c1, c2} = {10^6 C1, 10^(-11) C2} /. nlm["BestFitParameters"];
Plot[model[c1 , c2 , 100000, tmax][t], {t, 0, tmax},
Epilog :> Point[data], PlotRange -> All]


And you can make better predictions just fitting a line to a log-log plot:

mle = NonlinearModelFit[Log[data], a + b x, {a, b}, x]
Show[ListPlot[data],
Plot[Exp[mle[Log[x]]], {x, Min[data[[All, 1]]], Max[data[[All, 1]]]}]]


You'll need to decide if you want to fit the data or figure out why the data doesn't fit the theory.

• Thank you. It's OK that the data doesn't fit this theory. It shouldn't. I just wanted to eliminate mistakes before I start using more complicated theory with more differential equations. Also, that's why I didn't want to reformulate this simple PDE.
– K.J.
Sep 11, 2018 at 16:31
• So I think what you want to do is use the much, much speedier approach that @AlexTrounev demonstrated for the PDE and then maybe use the scaling that I outlined for the fit.
– JimB
Sep 11, 2018 at 18:27

In your model, you do not need to integrate into z, since the solution of the diffusion equation is practically independent of z. Also, do not use multiple coefficients in the model, it will be enough to use two that do not affect the solution, but affect the function that we will fit into the data.

data = {{8102., 0.00029}, {16294., 0.00034}, {24486.,
0.000379}, {32678., 0.000416}, {40870., 0.000437}, {49062.,
0.00045000000000000004}, {57254., 0.000471}, {66678.93400000001,
0.000486}, {87158.93400000001, 0.000512}}; tmax = First@Last@data;
H = 0.020; spInit[z_] := 0; dq = 10^5;
funLoad[t_] := If[t >= 10^-8, -dq, 0];
solSP = NDSolveValue[{ D[sp[z, t], {z, 2}] == D[sp[z, t], t],
sp[z, 0] == spInit[z], sp[0, t] == funLoad[t],
sp[H, t] == funLoad[t]}, sp, {z, 0, H}, {t, 0, tmax},
Method -> "BDF"];

solSett[t_?NumericQ, k1_?NumericQ, k2_?NumericQ] := -(H/k2)*
solSP[H/2, k1*t];
Plot[solSett[t, 6.7*10^-5, 4.2*10^6], {t, 0, tmax},
Epilog :> Point[data], PlotRange -> All] // AbsoluteTiming


• The fit is independent of the data? Love it. That's some theory!
– JimB
Sep 11, 2018 at 16:43
• @JimB What are you talking about? Sep 11, 2018 at 16:50
• If you set tmax = 87158.934, then what you show completely ignores the data. Although...I guess I assumed that you were providing a complete answer. Is the fitting part just left out or not there because it is unnecessary?
– JimB
Sep 11, 2018 at 17:01
• @JimB I answer the question: "Also I wonder: what can I do to speed up the calculations?". I showed that we can simplify the model, which works 488 times faster than the author's code. Sep 11, 2018 at 17:10
• Got it. Sorry for my confusion.
– JimB
Sep 11, 2018 at 17:16