# Plotting gradient and newton directions for parametric system of nonlinear ODEs

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;
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]}];