I have a bucket (system) of chemical kinetics models, a nonlinear dynamical system given by:
where kf >= 0 and kr >= 0 are the parameters. The initial conditions are A(0) = B(0) = 1 and C(0) = 0. I've generated data according to y1 = C(0.5)+noise and y2 = C(2)+noise where the noise is normally distributed with mu=0 and sigma=0.1 using kf = 0.1 and kr = 2.
odes = {A'[t] == -kf A[t] B[t] + kr C1[t],
B'[t] == -kf A[t] B[t] + kr C1[t],
C1'[t] == kf A[t] B[t] - kr C1[t], A[0] == 1, B[0] == 1,
C1[0] == 0};
odesData = odes /. {kf -> 0.1, kr -> 2};
soln = NDSolve[odesData, {A, B, C1}, {t, 0, 5}][[1]];
SeedRandom[10]
y1 = C1[0.5] + RandomVariate[NormalDistribution[0, 0.1]] /. soln
SeedRandom[30]
y2 = C1[2] + RandomVariate[NormalDistribution[0, 0.1]] /. soln
data = {{0.5, y1}, {2, y2}};
I need to visualize the data using an input/output picture, a parameter space picture, and a data space picture. For the parameter space and data space pictures, I also need to plot the gradient and newton directions for several randomly chosen parameter values.
Using ParametricNDSolve I am able to get an input/output picture.
kfmax = 1;
krmax = 5;
numsteps = 100;
kfrange = Range[0, kfmax, kfmax/numsteps];
krrange = Range[0, krmax, krmax/numsteps];
soln = ParametricNDSolve[odes, {A, B, C1}, {t, 0, 5}, {kf, kr}];
eqns = Evaluate[Table[C1[kf, kr][t] /. soln, {kf, kfrange}, {kr, krrange}]];
feweqns = Flatten[eqns][[ ;; ;; 1000]];
eqnsplot = Plot[feweqns, {t, 0, 2.5}, PlotRange -> All]
bestfitplot = Plot[model[kf, kr][t] /. fit, {t, 0, 2.5}, PlotRange -> All,
AxesLabel -> {"t", "C[t]"}];
Using ParametricNDSolveValue and FindFit I am able to get best fit parameters for the model (line is with fit parameters, and points are generated data).
model = ParametricNDSolveValue[odes, C1, {t, 0, 5}, {kf, kr}]
fit = FindFit[data, {model[kf, kr][t], {kf > 0, kr >= 0}}, {kf, kr}, t]
odes = {A'[t] == -kf A[t] B[t] + kr C1[t], B'[t] == -kf A[t] B[t] + kr C1[t],
C1'[t] == kf A[t] B[t] - kr C1[t], A[0] == 1, B[0] == 1, C1[0] == 0};
bestfitplot = Plot[model[kf, kr][t] /. fit, {t, 0, 2.5}, PlotRange -> All, AxesLabel -> {"t", "C[t]"}];
validptplot = ListPlot[{{0.5, y1}, {2, y2}}]; (*the data you're trying to fit*)
Show[bestfitplot, validptplot]
I am also able to visualize the parameter space using the log of the cost function:
list1 = eqns /. t -> 0.5;
list2 = eqns /. t -> 2;
cost = (list1 - y1)^2 + (list2 - y2)^2; cost // MatrixForm; (*kf-th row and kr-th column*)
parspaceplot = ListContourPlot[Log[cost], PlotLegends -> Automatic, DataRange -> {{0, krmax}, {0, kfmax}}, FrameLabel -> {"kr", "kf"}, Contours -> 50];
bestfitptplot = ListPlot[{{kr, kf}} /. fit];
parspace = Show[parspaceplot, bestfitptplot]
as well as the data space picture
max = 2;
dy1 = 0.2; dy2 = dy1;
modely1 = MapThread[model[#1, #2][0.5] &, Table[{i, j}, {i, 0, max, dy1}, {j, 0, max, dy2}][[#]]\[Transpose]] & /@ Range[max/dy1];
modely2 = MapThread[model[#1, #2][2] &, Table[{i, j}, {i, 0, max, dy1}, {j, 0, max, dy2}][[#]]\[Transpose]] & /@ Range[max/dy1];
plot1 = ListPlot[{modely1[[#]], modely2[[#]]}\[Transpose] & /@ Range[max/dy1], Joined -> True, PlotStyle -> Blue, AxesLabel -> {"y1", "y2"}];
plot2 = ListPlot[{modely1\[Transpose][[#]], modely2\[Transpose][[#]]}\[Transpose] & /@ Range[max/dy1], Joined -> True, PlotStyle -> Red];
plot3 = ListPlot[{{{y1, y2}}, {{model[kf, kr][0.5] /. fit, model[kf, kr][2] /. fit}}}, PlotStyle -> {{Black, Dot, PointSize[0.025]}, {Green}}];
Show[plot1, plot2, plot3]
However, when I try to compute the gradient (j'.r) and newton (solve (j'j).x==-grad for x) directions, where j === jacobian, j'j === fisher information matrix, r === residuals, and j' denotes the transpose of j, I don't get what I think I should be getting. I.e. the gradient direction isn't perpendicular to the contour lines.
SeedRandom[]
randmax = 10;
randoms = Table[{RandomReal[{0, krmax}], RandomReal[{0, kfmax}]}, {i,randmax}];
r = {y1, y2} - (model[kf, kr][#] & /@ {0.5, 2});
j = -(Grad[model[kf, kr][#], {kf, kr}] & /@ {0.5, 2})/0.1;
grad = j\[Transpose].r;
fish = j\[Transpose].j /. {kr -> randoms[[#, 1]], kf -> randoms[[#, 2]]} & /@ Range[randmax];
grads = grad /. {kr -> randoms[[#, 1]], kf -> randoms[[#, 2]]} & /@Range[randmax];
newts = LinearSolve[fish[[#]], -grads[[#]]] & /@ Range[randmax];
gradarrows = Graphics[{Black, Arrow[{randoms[[#]], randoms[[#]] +Normalize[grads[[#]]]/2.5}] & /@ Range[randmax]}];
newtarrows = Graphics[{Blue, Arrow[{randoms[[#]], randoms[[#]] + Normalize[newts[[#]]]/10}] & /@Range[randmax]}];
Show[parspace, gradarrows, newtarrows]
Things still don't look very good even if I do just cost instead of Log[cost].
This is a homework question for a graduate course in predictive modeling. A jupyter notebook is provided for the class and I know how I can numerically compute the jacobian (derive sensitivity equations and use odeint to solve the 9 equations simultaneously), but I got invested enough in this code that I wanted to see it to the end. I could do the same using Mathematica's NDSolve with the 9 equations, but it seemed like I should be able to get the jacobian using the parametric equation produced by ParametricNDSolve (similar to what is done in the parameter sensitivity topic in ParametricNDSolve's reference page).
Any suggestions on how I can get the gradient and newton directions? (Mathematical or coding advice both welcome).
P.s. This is my first post and the first time I wasn't able to resolve things with what I could find online and in this forum, which has been very helpful!
