See Figure 1 of
Zvi Hochman, Francois. Waldner
A simple N calculator for achieving water-limited yield of wheat crops
https://agronomyaustraliaproceedings.org/images/sampledata/2022/NutrientsAndFertility/ASAhochman_z_523s.pdf
The caption to their Figure 1 is as follows.
Wheat grain yield as a function of Available N (Nf) and Evapotranspiration (ET) for
(a) observed and (b) APSIM-simulated data.
Nf in kg/ha; Y in kg/ha; ET in mm
For the observed data yield was expressed as
Y = -1523 + 9 x Nf - 0.029 x Nf^2 + 16 x ET - 0.023 x ET^2 + 0.03 x ET x Nf;
(R^2 = 0.47; RMSE = 986 kg/ha; P < 0.001; Nf = 960).
For the APSIM simulated data yield was expressed as
Y = -1161 + 3 x Nf - 0.05 x Nf^2 + 14 x ET - 0.02 x ET^2 + 0.09 x ET x Nf;
(R^2 = 0.73; RMSE = 681 kg/ha; P < 0.001; Nf = 1,814).
(APSIM - from CSIRO and various State Ag. Depts - is one of several free open-source plant growth simulation packages.
There are firms that run APSIM or provide nicer online front-ends to it.
CropProphet is associated with Prescient Weather and seems as much involved with commodity trading as the actual farming.)
A bit about soil chemistry: soilchemn
Conjecture: yield is proportional to biomass (well some function of final biomass).
If we had a similar quadratic relation for biomass say at day x, say 20, (of the 180 day growth period), giving B[x] from Nf[0]
and them followed that up with another quadratic relation associated with application Nf[x],
perhaps (Y/c=) B[180]*quadratic in Nf[x]
one could optimise things like bf=B[180] over (Nf[0],Nf[x])
(Or perhaps priceB *B[180]*- costN*(Nf[0]+Nf[x]) )
Biquadratics occur a lot elsewhere. See biquadraticsgen.
That arising, as below, in the nitrogen-wheat application, seems in contour plots to look a bit like the twisted Edwards curve.
See wikipedia.
Edwards curve and related are well known amongst Comp. Sci students (about 3rd year undergrad level) because of their use in cryptography.
Also, it seems to GK an interesting challenge, given conditions on the signs of variables determining the coefficients of a biquadratic, to classify the critical points, give general info on them, how many and of what type.
Return to the nitrogen application. Lots of coefficients to determine, and this would involve growth of plant matters, notforGK!!
But after that it is very very simple. The lines starting check return 0.
cn is the ratio of costN to priceB
(To get the real profit, multiply the "profit" in the code by priceB.)
(* Mathematica, but sympy, numpy could do enough
January2026 version See elsewhere for later versions *)
bx = (a00 + a01*n0 + a02*n0^2);
bf = bx*(ax0 + ax1*nx + ax2*nx^2);
(* a00, a01,ax0,ax1 all positive; a02,ax2 both negative *)
profit = bf - cn*(n0 + nx);
(* profit is a biquadratic in (n0,nx) *)
(* The matrix items associated with the biquadratic are NOT used yet,
but are simple and pretty*)
M = CoefficientList[profit, {n0, nx}];
check = Simplify[({{1, n0, n0^2}} . M . Map[List, {1, nx, nx^2}])[[1, 1]] - profit]
Mcn = {{0, 1, 0}, {1, 0, 0}, {0, 0, 0}};
Mpbf = {{a00*ax0, a00*ax1, a00*ax2},
{a01*ax0, a01*ax1, a01*ax2},
{a02*ax0, a02*ax1, a02*ax2}};
check = Simplify[Mpbf - Mcn*cn - M]
(* Mpbf has a 2 dimensional nullspace and an eigenvalue lam *)
lam = a00*ax0 + a01*ax1 + a02*ax2;
evec ={ {a00,a01,a02}}
check = Simplify[Outer[Times, {a00, a01, a02}, {ax0, ax1, ax2}]- Mpbf]
(* vectors in the 2d nullspace
{-ax2, 0, ax0}, {-ax1, ax0, 0} and {0, -ax2, ax1} *)
(* NMaximize of profit with positivity constraints on n0 and nx
would solve the problem.
If we guess the optimum is not at boundary,both n0 >0 and nx>0,
can go for a calculus maximum, which we now do. *)
dprofit0 = Collect[D[profit, n0], n0, Factor]
dprofitx = Collect[D[profit, nx], nx, Factor]
hessprofit = D[profit, {{n0, nx}, 2}]
(* {{2 a02 (ax0 + ax1 nx + ax2 nx^2), (a01 + 2 a02 n0) (ax1 + 2 ax2 nx)},
{(a01 + 2 a02 n0) (ax1 + 2 ax2 nx), 2 ax2 (a00 + a01 n0 + a02 n0^2)}} *)
(* When cn=0, cost of N is 0, all the critical points, zeros of the gradient
are easily found. Two of them have both n0 and nx positive.
Obvious is that
{n0 -> -(a01/(2 a02)), nx -> -(ax1/(2 ax2))} is the local maximum
The other with both n0 and nx positive is a saddle. *)
(* If the cost of N is small, positive one uses a bit less *)
AsymptoticSolve[{dprofit0 == 0, dprofitx == 0} ,
{{n0, nx}, {-a01/(2*a02), -ax1/(2*ax2)}}, {cn, 0, 2}
(* {{n0 -> -(a01/(2 a02)) - (2 ax2 cn)/(a02 (ax1^2 - 4 ax0 ax2)),
nx -> -(ax1/(2 ax2)) - (2 a02 cn)/((a01^2 - 4 a00 a02) ax2)}} *)
(*General cn >0 :
Finding (n0,nx) jointly for max profit requires solving a quintic for n0 *)
Eliminate[{dprofit0 == 0, dprofitx == 0}, nx]
(*tidied a bit *)
n0Eq=cn^2 (a01 + 2 a02 n0) - 4 ax2 cn (a00 + a01 n0 + a02 n0^2)^2 -
(ax1^2 - 4 ax0 ax2) (a01 + 2 a02 n0) (a00 + a01 n0 + a02 n0^2)^2;
(* n0Eq has a quadratic factor when (ax0,ax1,ax2)=w*(a00,a01,a02)
The quadratic factor is
cn (a01 + 2 a02 n0) - (-a01^2 + 4 a00 a02) (a00 + a01 n0 + a02 n0^2)
*)
(* Once one has, numerically solved the quintic for n0, the equation
dprofitx == 0 gives nx *)
n0EqFn[n0_] := cn^2 (a01 + 2 a02 n0) -
4 ax2 cn (a00 + a01 n0 + a02 n0^2)^2 -
(ax1^2 - 4 ax0 ax2) (a01 + 2 a02 n0) (a00 + a01 n0 + a02 n0^2)^2;
Factor[n0EqFn[-(a01/(2 a02))]]
(* -(((-a01^2 + 4 a00 a02)^2 ax2 cn)/(4 a02^2)) nonnegative *)
FullSimplify[n0EqFn[(-a01 + Sqrt[a01^2 - 4 a00 a02])/(2 a02)]]
(* Sqrt[a01^2 - 4 a00 a02] cn^2 nonnegative *)
FullSimplify[n0EqFn[(-a01 - Sqrt[a01^2 - 4 a00 a02])/(2 a02)]]
(* -Sqrt[a01^2 - 4 a00 a02] cn^2 nonpositive *)
Factor[n0EqFn1[-(a01/(2 a02)) - (2 ax2 cn)/(a02 (ax1^2 - 4 ax0 ax2))]]
(* (4 ax2 cn^3)/(-ax1^2 + 4 ax0 ax2) nonnegative *)
Other things to look at
(even if biomass goes up with N, there are "diminishing returns").
Mewa Singh Dhanoa, Ruth Sanderson, Laura Maritza Cardenas
Available via ResearchGate
On the application of the Mitscherlich-Baule equation to fertiliser response data
November 2022CABI Reviews 2022(028)
DOI: 10.1079/cabireviews202217050