# Administrative Information *This is a supplement to the initial-value-problem unit (Day 15 onward). It records, at a careful graduate level, the complete multivariable Taylor-series derivation of the classical four-stage, fourth-order Runge–Kutta method for the scalar nonautonomous problem $y'=f(t,y)$, $y(t_n)=y_n$. The derivation of RK2 by the same technique is assumed known. Nothing here is quoted from rooted-tree theory; the eight order conditions are derived directly by expansion and coefficient matching, and the tree correspondence is noted only in a postscript. Every material algebraic step is shown.* > [!info]- Scope and prerequisites > - Assumed: the RK2 (midpoint/Heun family) derivation via second-order Taylor matching. > - Standing smoothness assumption: $f\in C^4$ on an open set containing the solution graph, so that all partial derivatives through total order four exist, are continuous, and may be evaluated freely; the exact solution is then $C^5$ in $t$. > - All expansions are in the step size $h>0$, about the base point $(t_n,y_n)$, which we abbreviate $(t,y)$. ## Notation All partial derivatives of $f$ are evaluated at the base point $(t,y)=(t_n,y_n)$ unless explicitly written otherwise. Subscripts denote partial differentiation: $ f_t=\frac{\partial f}{\partial t},\qquad f_y=\frac{\partial f}{\partial y},\qquad f_{tt}=\frac{\partial^2 f}{\partial t^2},\qquad f_{ty}=\frac{\partial^2 f}{\partial t\,\partial y},\qquad f_{yy}=\frac{\partial^2 f}{\partial y^2}, $ and at third order $ f_{ttt},\qquad f_{tty},\qquad f_{tyy},\qquad f_{yyy}, $ with the analogous meanings. Since $f\in C^4$, mixed partials are symmetric ($f_{ty}=f_{yt}$, etc.), and we always write the $t$-derivatives first. The unadorned symbol $f$ means $f(t_n,y_n)$. For a differentiable function $g(t,y)$ evaluated along the exact solution $y(t)$, the chain rule gives the *total derivative along the solution* $ \frac{d}{dt}\,g\bigl(t,y(t)\bigr)=g_t+g_y\,y'(t)=g_t+g_y f, $ the last equality holding because $y'=f$ along the solution. We will use this repeatedly, each time writing out the result rather than hiding it in operator notation. Four grouped quantities recur constantly, and it pays to name them now (each is a specific polynomial in the partials of $f$, evaluated at the base point): $ \begin{aligned} F &:= f_t+f\,f_y, \\ P &:= f_{ty}+f\,f_{yy}, \\ G &:= f_{tt}+2f\,f_{ty}+f^2 f_{yy}, \\ H &:= f_{ttt}+3f\,f_{tty}+3f^2 f_{tyy}+f^3 f_{yyy}. \end{aligned} $ $F$ is the total derivative of $f$ along the solution; $P$ is the total derivative of $f_y$ along the solution in the same sense ($P=(f_y)_t+(f_y)_y\,f$); $G$ and $H$ are the "pure" second and third directional derivatives of $f$ in the direction $(1,f)$: $ G=\Bigl(\partial_t+f\,\partial_y\Bigr)^{(2)}f \quad\text{and}\quad H=\Bigl(\partial_t+f\,\partial_y\Bigr)^{(3)}f $ *where the superscript means the directional derivative is taken with the direction vector $(1,f)$ held frozen at its base-point value* — that is, $G$ and $H$ are what you get from the binomial pattern $\sum_k \binom{n}{k} f^k\, \partial_t^{n-k}\partial_y^k f$, not from iterating the non-constant-coefficient operator $\partial_t + f(t,y)\partial_y$. This distinction matters: iterating the full operator would differentiate the coefficient $f$ as well, producing extra terms. Those extra terms are exactly the ones we will track separately below. --- # 1. The exact Taylor expansion of $y(t+h)$ through $O(h^4)$ Since $f\in C^4$, the solution $y(\cdot)$ is five times continuously differentiable, and Taylor's theorem with Lagrange remainder gives $ y(t+h)=y+h\,y'+\frac{h^2}{2}\,y''+\frac{h^3}{6}\,y'''+\frac{h^4}{24}\,y^{(4)}+\frac{h^5}{120}\,y^{(5)}(\xi), \qquad \xi\in(t,t+h). $ We must express $y',y'',y''',y^{(4)}$ in terms of partial derivatives of $f$ at $(t,y)$, by repeated total differentiation of $y'=f(t,y(t))$. ## 1.1 First derivative By the differential equation itself, $ y'=f. $ ## 1.2 Second derivative Differentiate $y'(t)=f(t,y(t))$ once in $t$ using the chain rule: $ y''=\frac{d}{dt}f(t,y(t))=f_t+f_y\,y'=f_t+f_y f=F. $ ## 1.3 Third derivative Differentiate $y''=f_t+f_y f$ term by term, applying $\frac{d}{dt}g=g_t+g_y f$ to each factor. First term: $ \frac{d}{dt}f_t=(f_t)_t+(f_t)_y\,f=f_{tt}+f_{ty}f. $ Second term, by the product rule: $ \frac{d}{dt}\bigl(f_y f\bigr) =\Bigl(\frac{d}{dt}f_y\Bigr)f+f_y\Bigl(\frac{d}{dt}f\Bigr) =\bigl(f_{ty}+f_{yy}f\bigr)f+f_y\bigl(f_t+f_y f\bigr). $ Expanding the products: $ \frac{d}{dt}\bigl(f_y f\bigr)=f\,f_{ty}+f^2 f_{yy}+f_y f_t+f\,f_y^2. $ Adding the two pieces and collecting: $ y'''=f_{tt}+2f\,f_{ty}+f^2 f_{yy}+f_y f_t+f\,f_y^2. $ In the grouped notation, since $G=f_{tt}+2ff_{ty}+f^2f_{yy}$ and $f_yF=f_yf_t+ff_y^2$, $ y'''=G+f_y F. $ ## 1.4 Fourth derivative Differentiate each of the five monomials in $y'''=f_{tt}+2f\,f_{ty}+f^2 f_{yy}+f_y f_t+f\,f_y^2$ separately. Throughout, $\frac{d}{dt}$ applied to any partial derivative $X$ of $f$ gives $X_t + X_y f$. **(i)** $\dfrac{d}{dt}f_{tt}=f_{ttt}+f_{tty}f.$ **(ii)** Product rule on $2f\,f_{ty}$: $ \frac{d}{dt}\bigl(2f\,f_{ty}\bigr) =2\bigl(f_t+f_yf\bigr)f_{ty}+2f\bigl(f_{tty}+f_{tyy}f\bigr) =2f_t f_{ty}+2f\,f_y f_{ty}+2f\,f_{tty}+2f^2 f_{tyy}. $ **(iii)** Product rule on $f^2 f_{yy}$, using $\frac{d}{dt}f^2 = 2f(f_t+f_yf)$: $ \frac{d}{dt}\bigl(f^2 f_{yy}\bigr) =2f\bigl(f_t+f_y f\bigr)f_{yy}+f^2\bigl(f_{tyy}+f_{yyy}f\bigr) =2f\,f_t f_{yy}+2f^2 f_y f_{yy}+f^2 f_{tyy}+f^3 f_{yyy}. $ **(iv)** Product rule on $f_y f_t$: $ \frac{d}{dt}\bigl(f_y f_t\bigr) =\bigl(f_{ty}+f_{yy}f\bigr)f_t+f_y\bigl(f_{tt}+f_{ty}f\bigr) =f_t f_{ty}+f\,f_t f_{yy}+f_y f_{tt}+f\,f_y f_{ty}. $ **(v)** Product rule on $f\,f_y^2$, using $\frac{d}{dt}f_y^2=2f_y(f_{ty}+f_{yy}f)$: $ \frac{d}{dt}\bigl(f\,f_y^2\bigr) =\bigl(f_t+f_yf\bigr)f_y^2+f\cdot 2f_y\bigl(f_{ty}+f_{yy}f\bigr) =f_y^2 f_t+f\,f_y^3+2f\,f_y f_{ty}+2f^2 f_y f_{yy}. $ Now sum (i)–(v) and collect like monomials. Listing every raw term: $ \begin{array}{lll} \text{from (i):} & f_{ttt}, & f\,f_{tty}\\[2pt] \text{from (ii):} & 2f_t f_{ty},\quad 2f\,f_y f_{ty}, & 2f\,f_{tty},\quad 2f^2 f_{tyy}\\[2pt] \text{from (iii):} & 2f\,f_t f_{yy},\quad 2f^2 f_y f_{yy}, & f^2 f_{tyy},\quad f^3 f_{yyy}\\[2pt] \text{from (iv):} & f_t f_{ty},\quad f\,f_t f_{yy}, & f_y f_{tt},\quad f\,f_y f_{ty}\\[2pt] \text{from (v):} & f_y^2 f_t,\quad f\,f_y^3, & 2f\,f_y f_{ty},\quad 2f^2 f_y f_{yy} \end{array} $ Collecting coefficients monomial by monomial: $ \begin{aligned} y^{(4)} ={}& f_{ttt} +3f\,f_{tty} +3f^2 f_{tyy} +f^3 f_{yyy} && (1+2,\;2+1 \text{ from (i)–(iii)})\\ &+f_y f_{tt} +5f\,f_y f_{ty} +4f^2 f_y f_{yy} && (2+1+2 \text{ and } 2+2 \text{ from (ii),(iv),(v)})\\ &+3f_t f_{ty} +3f\,f_t f_{yy} && (2+1 \text{ from (ii),(iii),(iv)})\\ &+f_y^2 f_t +f\,f_y^3. && (\text{from (v)}) \end{aligned} $ This is eleven distinct monomials. We now verify the grouped form. Compute each group: $ f_y G=f_y f_{tt}+2f\,f_y f_{ty}+f^2 f_y f_{yy}, $ $ 3FP=3\bigl(f_t+f f_y\bigr)\bigl(f_{ty}+f f_{yy}\bigr) =3f_t f_{ty}+3f\,f_t f_{yy}+3f\,f_y f_{ty}+3f^2 f_y f_{yy}, $ $ f_y^2 F=f_y^2 f_t+f\,f_y^3. $ Adding $H+f_yG+3FP+f_y^2F$ and collecting: the $f\,f_yf_{ty}$ coefficient is $2+3=5$; the $f^2f_yf_{yy}$ coefficient is $1+3=4$; every other monomial appears in exactly one group with the coefficient listed above. The totals agree monomial-for-monomial with the sum (i)–(v), so $ \boxed{\,y^{(4)}=H+f_y G+3FP+f_y^2 F\,} $ ## 1.5 The exact expansion Substituting $y'=f$, $y''=F$, $y'''=G+f_yF$, $y^{(4)}=H+f_yG+3FP+f_y^2F$ into the Taylor polynomial: $ y(t+h)=y+hf+\frac{h^2}{2}F+\frac{h^3}{6}\bigl(G+f_yF\bigr) +\frac{h^4}{24}\bigl(H+f_yG+3FP+f_y^2F\bigr)+O(h^5). $ **What has been established.** The exact solution admits the expansion above, with every coefficient computed by explicit chain- and product-rule differentiation; the $O(h^5)$ remainder is uniform on compact sets by continuity of the fifth derivative (equivalently, of the fourth partials of $f$). The eleven fourth-order monomials and their multiplicities $1,3,3,1,1,5,4,3,3,1,1$ have been derived, not quoted. --- # 2. The elementary differentials through order four By an *elementary differential of order $q$* we mean one of the specific polynomial expressions in $f$ and its partials that can appear in $y^{(q)}$; these are the building blocks against which both expansions will be matched. From §1: **Order 1** (1 term): $ f. $ **Order 2** (2 monomials, one group): $ f_t,\qquad f\,f_y \qquad\text{grouped as } F=f_t+f f_y. $ **Order 3** (5 monomials, two groups): $ f_{tt},\quad f\,f_{ty},\quad f^2 f_{yy} \qquad\text{grouped as } G, $ $ f_y f_t,\quad f\,f_y^2 \qquad\text{grouped as } f_y F. $ **Order 4** (11 monomials, four groups): $ \begin{array}{ll} \text{group } H: & f_{ttt},\quad f\,f_{tty},\quad f^2 f_{tyy},\quad f^3 f_{yyy}\\[3pt] \text{group } f_yG: & f_y f_{tt},\quad f\,f_y f_{ty},\quad f^2 f_y f_{yy}\\[3pt] \text{group } FP: & f_t f_{ty},\quad f\,f_t f_{yy},\quad f\,f_y f_{ty},\quad f^2 f_y f_{yy}\\[3pt] \text{group } f_y^2F: & f_y^2 f_t,\quad f\,f_y^3 \end{array} $ Two cautions, both of which matter later and justify the instruction not to combine terms prematurely. *First*, the groups overlap at the monomial level: $f\,f_yf_{ty}$ occurs in both $f_yG$ (with weight $2$) and $FP$ (with weight $1$), and $f^2f_yf_{yy}$ occurs in both $f_yG$ (weight $1$) and $FP$ (weight $1$). So "matching coefficients of $f_yG$ and of $FPquot; is meaningful only if we can first show the four groups are independent — which we do by exhibiting *marker monomials*: $f_{ttt}$ occurs only in $H$; $f_yf_{tt}$ only in $f_yG$; $f_tf_{ty}$ only in $FP$; $f\,f_y^3$ only in $f_y^2F$. Matching will therefore proceed at the monomial level, with the markers pinning the group coefficients and the shared monomials serving as consistency checks (§7). *Second*, the groups — not the monomials — are the natural objects, because (as we will see) every stage expansion produces the partials of $f$ only in the combinations $F,P,G,H$ once the internal-consistency relations hold. The four order-4 groups are exactly the four elementary differentials of rooted-tree theory specialized to the scalar nonautonomous case, but we make no use of that fact until the postscript. **What has been established.** A finite list of $1+2+5+11=19$ monomials spans everything that can appear through $h^4$ in either expansion, organized into $1+1+2+4=8$ groups whose independence is certified by marker monomials. --- # 3. The stage expansions: general principles, $k_1$, $k_2$, and the origin of the consistency relations ## 3.1 The ansatz and the accuracy budget The general explicit four-stage method is $ \begin{aligned} k_1 &= f(t,y),\\ k_2 &= f\bigl(t+c_2h,\;y+h\,a_{21}k_1\bigr),\\ k_3 &= f\bigl(t+c_3h,\;y+h(a_{31}k_1+a_{32}k_2)\bigr),\\ k_4 &= f\bigl(t+c_4h,\;y+h(a_{41}k_1+a_{42}k_2+a_{43}k_3)\bigr), \end{aligned} \qquad y_{n+1}=y+h\bigl(b_1k_1+b_2k_2+b_3k_3+b_4k_4\bigr). $ **Why each stage needs only $O(h^3)$ accuracy.** Every $k_i$ enters $y_{n+1}$ multiplied by $hb_i$. Our target is to match the exact expansion through $h^4$, i.e. to make the one-step (local truncation) error $O(h^5)$. If $k_i$ is known with an error of size $O(h^4)$, that error contributes $hb_i\cdot O(h^4)=O(h^5)$ to $y_{n+1}$ — invisible at the order we are matching. Hence it suffices to expand each stage exactly through the $h^3$ terms and carry a single symbol $O(h^4)$ for the rest. (Conversely, an error of size $O(h^3)$ in a stage would contaminate the $h^4$ terms of $y_{n+1}$, so nothing coarser suffices.) **The tool: two-variable Taylor expansion.** For offsets $\delta$ and $\Delta$ with $\delta=O(h)$, $\Delta=O(h)$, Taylor's theorem in two variables for $f\in C^4$ gives $ \begin{aligned} f(t+\delta,\,y+\Delta) ={}&f+\delta f_t+\Delta f_y\\ &+\tfrac12\bigl(\delta^2 f_{tt}+2\delta\Delta\,f_{ty}+\Delta^2 f_{yy}\bigr)\\ &+\tfrac16\bigl(\delta^3 f_{ttt}+3\delta^2\Delta\,f_{tty}+3\delta\Delta^2 f_{tyy}+\Delta^3 f_{yyy}\bigr) +R_4, \end{aligned} $ where the remainder $R_4$ collects the fourth-order terms with partials evaluated at an intermediate point; since $f\in C^4$ these are bounded, so $R_4=O\bigl((|\delta|+|\Delta|)^4\bigr)=O(h^4)$. All partials on the right are at the base point $(t,y)$. Every stage expansion below is an instance of this formula with $\delta=c_ih$ and $\Delta=$ the $y$-increment of that stage. ## 3.2 Stage one There is nothing to expand: $ k_1=f. $ ## 3.3 Stage two Here $\delta=c_2h$ and $\Delta=h\,a_{21}k_1=h\,a_{21}f$. Substitute into the two-variable formula, term by term: - Linear terms: $\;\delta f_t+\Delta f_y=h\bigl(c_2f_t+a_{21}f\,f_y\bigr)$. - Quadratic terms: $\;\tfrac12\bigl(c_2^2h^2f_{tt}+2c_2h\cdot a_{21}hf\cdot f_{ty}+a_{21}^2h^2f^2f_{yy}\bigr)=\tfrac{h^2}{2}\bigl(c_2^2f_{tt}+2c_2a_{21}f\,f_{ty}+a_{21}^2f^2f_{yy}\bigr)$. - Cubic terms: $\;\tfrac{h^3}{6}\bigl(c_2^3f_{ttt}+3c_2^2a_{21}f\,f_{tty}+3c_2a_{21}^2f^2f_{tyy}+a_{21}^3f^3f_{yyy}\bigr)$. So, exactly, $ k_2=f+h\bigl(c_2f_t+a_{21}f f_y\bigr) +\frac{h^2}{2}\bigl(c_2^2f_{tt}+2c_2a_{21}f f_{ty}+a_{21}^2f^2f_{yy}\bigr) +\frac{h^3}{6}\bigl(c_2^3f_{ttt}+3c_2^2a_{21}f f_{tty}+3c_2a_{21}^2f^2f_{tyy}+a_{21}^3f^3f_{yyy}\bigr)+O(h^4). $ ## 3.4 Where the internal-consistency relations come from Notice what the $h^1$ term of $k_2$ looks like: $c_2f_t+a_{21}f f_y$. The exact solution's second derivative is $F=f_t+ff_y$ — the two monomials $f_t$ and $ff_y$ occur there *with equal weight*. In the Runge–Kutta expansion, however, the time offset $c_2h$ feeds the $t$-partials and the state offset $h a_{21}f$ feeds the $y$-partials, and these are controlled by *different* free parameters. The same decoupling recurs at every order: $c_i$ powers multiply the $t$-heavy monomials and row sums of $A$ multiply the $y$-heavy ones. Consequence: if we match coefficients monomial-by-monomial with $c_2\neq a_{21}$ left independent, then already at order $h^2$ of $y_{n+1}$ we get *two* conditions, $ \sum_i b_ic_i=\tfrac12 \quad(\text{from } f_t), \qquad \sum_i b_i\Bigl(\textstyle\sum_j a_{ij}\Bigr)=\tfrac12 \quad(\text{from } f f_y), $ and at order $h^3$ the group $G$ splits into *three* conditions (from $f_{tt}$, $f f_{ty}$, $f^2f_{yy}$ respectively): $ \sum_i b_ic_i^2=\tfrac13,\qquad \sum_i b_ic_i\Bigl(\textstyle\sum_j a_{ij}\Bigr)=\tfrac13,\qquad \sum_i b_i\Bigl(\textstyle\sum_j a_{ij}\Bigr)^2=\tfrac13, $ and so on with worsening multiplicity at order four. All of these collapse pairwise (and triple-wise) into single conditions precisely when $ c_i=\sum_j a_{ij} \qquad\Longleftrightarrow\qquad c_2=a_{21},\quad c_3=a_{31}+a_{32},\quad c_4=a_{41}+a_{42}+a_{43}, $ because then every mixed power $c_i^{p}\bigl(\sum_ja_{ij}\bigr)^{q}$ equals $c_i^{p+q}$. The relations have an equivalent geometric reading. The $i$-th stage evaluates $f$ at the point $\bigl(t+c_ih,\;Y_i\bigr)$ with $Y_i=y+h\sum_ja_{ij}k_j$. Since every $k_j=f+O(h)$, $ Y_i=y+h\Bigl(\sum_j a_{ij}\Bigr)f+O(h^2), \qquad\text{while}\qquad y(t+c_ih)=y+c_ih\,f+O(h^2). $ So the sample point $(t+c_ih,Y_i)$ lies within $O(h^2)$ of the exact solution graph — i.e. the stage samples the slope field *where the solution actually is*, to leading order — if and only if $\sum_ja_{ij}=c_i$. Without the relation, each stage samples the slope field at a point displaced $O(h)$ horizontally relative to its vertical position on the graph, and the $t$- and $y$-expansions must then be matched separately, which is exactly the proliferation of conditions seen above. (One can also phrase this via the augmented autonomous system obtained by adjoining $t'=1$: applying the same tableau to the augmented system advances the $t$-component of stage $i$ by $h\sum_ja_{ij}\cdot 1$, while the nonautonomous format above uses $c_ih$; the two formats define the same method for all $f$ iff $c_i=\sum_ja_{ij}$. This is a repackaging of the previous paragraph, not an extra assumption: the augmented expansion is the same two-variable expansion with $t$ treated as a solution component.) We henceforth impose the consistency relations. We do not lose the classical method by doing so — its tableau satisfies them — and every subsequent expansion simplifies decisively. It must be conceded that this is a (standard, and universally adopted) *structural choice*: methods violating $c_i=\sum_ja_{ij}$ are not thereby excluded from attaining order four for the autonomous problem, but for the nonautonomous problem they would have to satisfy the larger split system, of which the collapsed system below is the tidy special case. ## 3.5 Stage two, simplified With $a_{21}=c_2$, the expansion of §3.3 becomes, term by term: - $h$: $\;c_2f_t+c_2ff_y=c_2\bigl(f_t+ff_y\bigr)=c_2F$. - $h^2$: $\;\tfrac12\bigl(c_2^2f_{tt}+2c_2^2ff_{ty}+c_2^2f^2f_{yy}\bigr)=\tfrac{c_2^2}{2}\bigl(f_{tt}+2ff_{ty}+f^2f_{yy}\bigr)=\tfrac{c_2^2}{2}G$. - $h^3$: $\;\tfrac{c_2^3}{6}\bigl(f_{ttt}+3ff_{tty}+3f^2f_{tyy}+f^3f_{yyy}\bigr)=\tfrac{c_2^3}{6}H$. $ \boxed{\;k_2=f+c_2h\,F+\frac{c_2^2h^2}{2}\,G+\frac{c_2^3h^3}{6}\,H+O(h^4).\;} $ Observe what is *absent*: $k_2$ contains no $f_yF$, no $f_yG$, no $FP$, no $f_y^2F$ terms. Those groups all involve $f_y$ acting on a *lower-order elementary differential*, and they can only arise when a stage increment itself contains $h^2$ or $h^3$ corrections built from earlier stages — that is, from *nesting*. $k_2s increment $hc_2f$ is proportional to the exact first-order increment, so $k_2$ is nothing but the beginning of the Taylor series of $f$ along the ray through $(t,y)$ with direction $(1,f)$ — which is precisely why only $F,G,H$ appear. **What has been established.** Each stage need only be expanded through $h^3$. The two-variable Taylor formula has been fixed once and for all. The consistency relations $c_i=\sum_ja_{ij}$ have been *derived* as the condition under which (a) the split $t$-versus-$y$ matching conditions collapse, equivalently (b) each stage samples the slope field within $O(h^2)$ of the solution graph; they are imposed from here on. And $k_2=f+c_2hF+\tfrac{c_2^2h^2}{2}G+\tfrac{c_2^3h^3}{6}H+O(h^4)$, with the nested groups provably absent. --- # 4. Full expansion of $k_3$ through $O(h^3)$ ## 4.1 The increment Stage three evaluates $f$ at $\bigl(t+c_3h,\;y+\Delta_3\bigr)$ with $ \Delta_3=h\bigl(a_{31}k_1+a_{32}k_2\bigr). $ Substitute $k_1=f$ and the boxed expansion of $k_2$. Since $\Delta_3$ is multiplied by at most $f_y$ (one power of an $h$-free partial) before entering an $h^3$ term of $k_3$, we need $\Delta_3$ itself through $h^3$: $ \Delta_3 =h\Bigl(a_{31}f+a_{32}\Bigl[f+c_2hF+\frac{c_2^2h^2}{2}G+O(h^3)\Bigr]\Bigr) =h\bigl(a_{31}+a_{32}\bigr)f+h^2a_{32}c_2F+h^3\,\frac{a_{32}c_2^2}{2}G+O(h^4). $ Using the consistency relation $a_{31}+a_{32}=c_3$: $ \Delta_3=hc_3f+h^2a_{32}c_2F+h^3\,\frac{a_{32}c_2^2}{2}G+O(h^4). $ For the quadratic and cubic Taylor terms we also need the powers of $\Delta_3$, again only through $h^3$: $ \Delta_3^2 =\Bigl(hc_3f+h^2a_{32}c_2F+O(h^3)\Bigr)^2 =h^2c_3^2f^2+2h^3c_3a_{32}c_2\,fF+O(h^4), $ where the cross term is $2\cdot(hc_3f)\cdot(h^2a_{32}c_2F)$ and the square of the $h^2$ part is $O(h^4)$; $ \Delta_3^3=\Bigl(hc_3f+O(h^2)\Bigr)^3=h^3c_3^3f^3+O(h^4). $ Finally the mixed products with $\delta_3=c_3h$: $ \delta_3\Delta_3=h^2c_3^2f+h^3c_3a_{32}c_2F+O(h^4), \qquad \delta_3^2\Delta_3=h^3c_3^3f+O(h^4), \qquad \delta_3\Delta_3^2=h^3c_3^3f^2+O(h^4). $ ## 4.2 Substitution into the two-variable expansion Insert $\delta_3=c_3h$ and the above into the master formula of §3.1, keeping everything through $h^3$. **Constant and linear part:** $ f+\delta_3f_t+\Delta_3f_y =f+hc_3f_t +\Bigl(hc_3f+h^2a_{32}c_2F+h^3\frac{a_{32}c_2^2}{2}G\Bigr)f_y+O(h^4), $ which regroups as $ f+hc_3\bigl(f_t+ff_y\bigr)+h^2a_{32}c_2\,f_yF+h^3\frac{a_{32}c_2^2}{2}\,f_yG+O(h^4) =f+hc_3F+h^2a_{32}c_2\,f_yF+h^3\frac{a_{32}c_2^2}{2}\,f_yG+O(h^4). $ **Quadratic part:** $\tfrac12\bigl(\delta_3^2f_{tt}+2\delta_3\Delta_3f_{ty}+\Delta_3^2f_{yy}\bigr)$ equals $ \tfrac12\Bigl(h^2c_3^2f_{tt} +2\bigl[h^2c_3^2f+h^3c_3a_{32}c_2F\bigr]f_{ty} +\bigl[h^2c_3^2f^2+2h^3c_3a_{32}c_2fF\bigr]f_{yy}\Bigr)+O(h^4). $ Separate the $h^2$ and $h^3$ contents: $ h^2\colon\quad \tfrac{c_3^2}{2}\bigl(f_{tt}+2ff_{ty}+f^2f_{yy}\bigr)=\tfrac{c_3^2}{2}G, $ $ h^3\colon\quad \tfrac12\bigl(2c_3a_{32}c_2F\,f_{ty}+2c_3a_{32}c_2\,fF\,f_{yy}\bigr) =c_3a_{32}c_2\,F\bigl(f_{ty}+ff_{yy}\bigr) =c_3a_{32}c_2\,FP. $ **Cubic part:** $\tfrac16\bigl(\delta_3^3f_{ttt}+3\delta_3^2\Delta_3f_{tty}+3\delta_3\Delta_3^2f_{tyy}+\Delta_3^3f_{yyy}\bigr)$; every factor already carries $h^3$, so only the leading terms of the powers survive: $ \tfrac{h^3}{6}\bigl(c_3^3f_{ttt}+3c_3^3ff_{tty}+3c_3^3f^2f_{tyy}+c_3^3f^3f_{yyy}\bigr)+O(h^4) =\frac{c_3^3h^3}{6}\,H+O(h^4). $ ## 4.3 Collected result $ \boxed{\; k_3=f+c_3h\,F +h^2\Bigl(\frac{c_3^2}{2}\,G+a_{32}c_2\,f_yF\Bigr) +h^3\Bigl(\frac{c_3^3}{6}\,H+c_3a_{32}c_2\,FP+\frac{a_{32}c_2^2}{2}\,f_yG\Bigr) +O(h^4).\;} $ Note the structural bookkeeping, which repeats at stage four: - The "ray" groups $F,G,H$ appear with the same coefficients $c_3,\tfrac{c_3^2}{2},\tfrac{c_3^3}{6}$ as in $k_2$ (with $c_3$ in place of $c_2$): these come from the leading part $hc_3f$ of the increment. - $f_yF$ at $h^2$ carries $a_{32}c_2$: the $h^2$ *correction* of the increment ($a_{32}$ weighting $k_2s $c_2hF$ term) hit by $f_y$. - $FP$ at $h^3$ carries $c_3\cdot a_{32}c_2$: the same $h^2$ correction interacting with the leading increment through the *second*-order partials. - $f_yG$ at $h^3$ carries $\tfrac{a_{32}c_2^2}{2}$: the $h^3$ correction of the increment (from $k_2s $G$ term) hit by $f_y$. - There is **no** $f_y^2F$ term. That group requires $f_y$ to act on an increment correction which *itself* contains an $f_yF$ term — i.e. double nesting. But $k_2$, feeding $\Delta_3$, has no $f_yF$ term. Double nesting first becomes possible at stage four, fed by $k_3$. **What has been established.** The complete $O(h^3)$ expansion of $k_3$, with each coefficient traced to its source in the increment algebra, and a structural explanation of which elementary-differential groups can and cannot appear. --- # 5. Full expansion of $k_4$ through $O(h^3)$ ## 5.1 The increment $ \Delta_4=h\bigl(a_{41}k_1+a_{42}k_2+a_{43}k_3\bigr). $ Substitute $k_1=f$, the boxed $k_2$ (§3.5), and the boxed $k_3$ (§4.3), each through $h^2$ (a further factor $h$ is supplied by the prefactor, and $\Delta_4$ is needed through $h^3$): $ \Delta_4=h\,a_{41}f +h\,a_{42}\Bigl[f+c_2hF+\frac{c_2^2h^2}{2}G\Bigr] +h\,a_{43}\Bigl[f+c_3hF+h^2\Bigl(\frac{c_3^2}{2}G+a_{32}c_2f_yF\Bigr)\Bigr]+O(h^4). $ Collect by powers of $h$: $ \Delta_4 =h\bigl(a_{41}+a_{42}+a_{43}\bigr)f +h^2\bigl(a_{42}c_2+a_{43}c_3\bigr)F +h^3\Bigl(\frac{a_{42}c_2^2+a_{43}c_3^2}{2}\,G+a_{43}a_{32}c_2\,f_yF\Bigr)+O(h^4). $ With the consistency relation $a_{41}+a_{42}+a_{43}=c_4$, and abbreviating $ \sigma:=a_{42}c_2+a_{43}c_3, \qquad \rho:=a_{42}c_2^2+a_{43}c_3^2, $ we have $ \Delta_4=hc_4f+h^2\sigma F+h^3\Bigl(\frac{\rho}{2}\,G+a_{43}a_{32}c_2\,f_yF\Bigr)+O(h^4). $ The needed powers, through $h^3$ (identical algebra to §4.1 with $c_4,\sigma$ in place of $c_3,\;a_{32}c_2$): $ \Delta_4^2=h^2c_4^2f^2+2h^3c_4\sigma\,fF+O(h^4), \qquad \Delta_4^3=h^3c_4^3f^3+O(h^4), $ $ \delta_4\Delta_4=h^2c_4^2f+h^3c_4\sigma F+O(h^4), \qquad \delta_4^2\Delta_4=h^3c_4^3f+O(h^4), \qquad \delta_4\Delta_4^2=h^3c_4^3f^2+O(h^4), $ where $\delta_4=c_4h$. Note that the $h^3$ part of $\Delta_4$ never enters the quadratic or cubic Taylor terms at this order — it survives only through the linear term $\Delta_4f_y$. ## 5.2 Substitution into the two-variable expansion **Constant and linear part:** $ f+\delta_4f_t+\Delta_4f_y =f+hc_4F +h^2\sigma\,f_yF +h^3\Bigl(\frac{\rho}{2}\,f_yG+a_{43}a_{32}c_2\,f_y^2F\Bigr)+O(h^4), $ where, exactly as in §4.2, the $h^1$ pieces combined via $c_4f_t+c_4ff_y=c_4F$, and the $h^3$ piece of $\Delta_4$ was multiplied by $f_y$ — this is the birth of the $f_y^2F$ group: $f_y$ acting on $\Delta_4s own $f_yF$ content, which $\Delta_4$ inherited from $k_3$, which in turn inherited it from $k_2s $F$ term. Three generations of stages are needed, which is why no three-stage explicit method can be fourth order (that claim is completed by the count of conditions in §9). **Quadratic part:** $ \tfrac12\Bigl(h^2c_4^2f_{tt} +2\bigl[h^2c_4^2f+h^3c_4\sigma F\bigr]f_{ty} +\bigl[h^2c_4^2f^2+2h^3c_4\sigma fF\bigr]f_{yy}\Bigr)+O(h^4), $ which splits as $ h^2\colon\;\frac{c_4^2}{2}\,G, \qquad h^3\colon\;\tfrac12\bigl(2c_4\sigma Ff_{ty}+2c_4\sigma fFf_{yy}\bigr)=c_4\sigma\,FP. $ **Cubic part:** $ \frac{h^3}{6}\bigl(c_4^3f_{ttt}+3c_4^3ff_{tty}+3c_4^3f^2f_{tyy}+c_4^3f^3f_{yyy}\bigr)+O(h^4) =\frac{c_4^3h^3}{6}\,H+O(h^4). $ ## 5.3 Collected result $ \boxed{\; k_4=f+c_4h\,F +h^2\Bigl(\frac{c_4^2}{2}\,G+\sigma\,f_yF\Bigr) +h^3\Bigl(\frac{c_4^3}{6}\,H+c_4\sigma\,FP+\frac{\rho}{2}\,f_yG+a_{43}a_{32}c_2\,f_y^2F\Bigr) +O(h^4),\;} $ with $\sigma=a_{42}c_2+a_{43}c_3$ and $\rho=a_{42}c_2^2+a_{43}c_3^2$. **What has been established.** The complete $O(h^3)$ expansions of all four stages. In generalized form, writing $d_i:=\sum_ja_{ij}c_j$, $e_i:=\sum_ja_{ij}c_j^2$, $g_i:=\sum_j\sum_l a_{ij}a_{jl}c_l$ (empty sums vanish for an explicit method), every stage obeys the single pattern $ k_i=f+c_ih\,F+h^2\Bigl(\frac{c_i^2}{2}G+d_i\,f_yF\Bigr) +h^3\Bigl(\frac{c_i^3}{6}H+c_id_i\,FP+\frac{e_i}{2}\,f_yG+g_i\,f_y^2F\Bigr)+O(h^4), $ with $(d_1,e_1,g_1)=(0,0,0)$, $(d_2,e_2,g_2)=(0,0,0)$, $(d_3,e_3,g_3)=(a_{32}c_2,\,a_{32}c_2^2,\,0)$, $(d_4,e_4,g_4)=(\sigma,\,\rho,\,a_{43}a_{32}c_2)$ — as direct inspection of the boxed formulas confirms. --- # 6. Assembling $y_{n+1}=y+h\sum_i b_ik_i$ Insert the four boxed stage expansions into the update and collect by powers of $h$. Since each $k_i$ is multiplied by $h$, the $h^0$ part of $k_i$ feeds the $h^1$ term of $y_{n+1}$, and so on. **Order $h^1$** (from the $f$ terms of all stages): $ h\bigl(b_1+b_2+b_3+b_4\bigr)f. $ **Order $h^2$** (from the $c_ihF$ terms; $k_1$ contributes nothing since $c_1=0$): $ h^2\bigl(b_2c_2+b_3c_3+b_4c_4\bigr)F. $ **Order $h^3$** (from the $h^2$ terms of the stages). The $G$ content: $ h^3\cdot\frac{1}{2}\bigl(b_2c_2^2+b_3c_3^2+b_4c_4^2\bigr)G, $ and the $f_yF$ content (only $k_3$ and $k_4$ possess one): $ h^3\bigl(b_3a_{32}c_2+b_4(a_{42}c_2+a_{43}c_3)\bigr)f_yF =h^3\Bigl(\sum_{i,j}b_ia_{ij}c_j\Bigr)f_yF. $ **Order $h^4$** (from the $h^3$ terms of the stages). Reading the four groups off the boxed formulas: $ H:\quad h^4\cdot\frac{1}{6}\bigl(b_2c_2^3+b_3c_3^3+b_4c_4^3\bigr), $ $ FP:\quad h^4\bigl(b_3c_3\,a_{32}c_2+b_4c_4\,(a_{42}c_2+a_{43}c_3)\bigr) =h^4\sum_{i,j}b_ic_ia_{ij}c_j, $ $ f_yG:\quad h^4\cdot\frac{1}{2}\bigl(b_3a_{32}c_2^2+b_4(a_{42}c_2^2+a_{43}c_3^2)\bigr) =h^4\cdot\frac12\sum_{i,j}b_ia_{ij}c_j^2, $ $ f_y^2F:\quad h^4\,b_4a_{43}a_{32}c_2 =h^4\sum_{i,j,l}b_ia_{ij}a_{jl}c_l. $ Altogether: $ \begin{aligned} y_{n+1}=y &+h\Bigl(\sum_ib_i\Bigr)f +h^2\Bigl(\sum_ib_ic_i\Bigr)F\\ &+h^3\biggl[\frac12\Bigl(\sum_ib_ic_i^2\Bigr)G+\Bigl(\sum_{i,j}b_ia_{ij}c_j\Bigr)f_yF\biggr]\\ &+h^4\biggl[\frac16\Bigl(\sum_ib_ic_i^3\Bigr)H +\Bigl(\sum_{i,j}b_ic_ia_{ij}c_j\Bigr)FP +\frac12\Bigl(\sum_{i,j}b_ia_{ij}c_j^2\Bigr)f_yG +\Bigl(\sum_{i,j,l}b_ia_{ij}a_{jl}c_l\Bigr)f_y^2F\biggr] +O(h^5). \end{aligned} $ **What has been established.** The full method expansion through $h^4$, organized by the same eight elementary-differential groups as the exact expansion of §1.5. Every coefficient is an explicit polynomial in the tableau entries, derived — not asserted — from the stage algebra of §§3–5. --- # 7. Matching coefficients: the eight order conditions Set the expansion of §6 equal, group by group, to the exact expansion $ y(t+h)=y+hf+\frac{h^2}{2}F+\frac{h^3}{6}G+\frac{h^3}{6}f_yF +\frac{h^4}{24}H+\frac{3h^4}{24}FP+\frac{h^4}{24}f_yG+\frac{h^4}{24}f_y^2F+O(h^5). $ **Why group-by-group matching is legitimate.** The two expansions are polynomial identities in the (freely choosable) values of $f$ and its partials at one point. Matching must therefore hold monomial by monomial across the 19 monomials of §2. At orders 1–3 each monomial belongs to exactly one group and (by §§3–5) appears inside that group with the same internal weights on both sides, so monomial matching is equivalent to group matching. At order 4, the marker monomials $f_{ttt}$, $f_yf_{tt}$, $f_tf_{ty}$, $ff_y^3$ each live in a single group and force the four group-coefficient equations below. The two shared monomials then need checking: on the method side the coefficient of $f\,f_yf_{ty}$ is $2\times(\text{coeff of }f_yG)+1\times(\text{coeff of }FP)$ and on the exact side it is $2\cdot\frac1{24}+1\cdot\frac{3}{24}=\frac{5}{24}$, matching §1.4's multiplicity $5$ of $y^{(4)}/24$; similarly $f^2f_yf_{yy}$ carries $1\cdot\frac1{24}+1\cdot\frac{3}{24}=\frac4{24}$ against multiplicity $4$. So once the four group equations hold, the shared monomials are automatically matched: no additional conditions hide in the overlaps. The matching, condition by condition, each labeled with its generating elementary differential: **Order 1, from $f$:** $ \text{(O1)}\qquad b_1+b_2+b_3+b_4=1. $ **Order 2, from $F=f_t+ff_y$:** $ \text{(O2)}\qquad b_2c_2+b_3c_3+b_4c_4=\frac12. $ **Order 3, from $G$** (method: $\tfrac12\sum b_ic_i^2$; exact: $\tfrac16$): $ \text{(O3)}\qquad b_2c_2^2+b_3c_3^2+b_4c_4^2=\frac13. $ **Order 3, from $f_yF$** (method: $\sum b_ia_{ij}c_j$; exact: $\tfrac16$): $ \text{(O4)}\qquad b_3a_{32}c_2+b_4\bigl(a_{42}c_2+a_{43}c_3\bigr)=\frac16. $ **Order 4, from $H$** (method: $\tfrac16\sum b_ic_i^3$; exact: $\tfrac1{24}$): $ \text{(O5)}\qquad b_2c_2^3+b_3c_3^3+b_4c_4^3=\frac14. $ **Order 4, from $FP$** (method: $\sum b_ic_ia_{ij}c_j$; exact: $\tfrac{3}{24}=\tfrac18$): $ \text{(O6)}\qquad b_3c_3a_{32}c_2+b_4c_4\bigl(a_{42}c_2+a_{43}c_3\bigr)=\frac18. $ **Order 4, from $f_yG$** (method: $\tfrac12\sum b_ia_{ij}c_j^2$; exact: $\tfrac1{24}$): $ \text{(O7)}\qquad b_3a_{32}c_2^2+b_4\bigl(a_{42}c_2^2+a_{43}c_3^2\bigr)=\frac1{12}. $ **Order 4, from $f_y^2F$** (method: $\sum b_ia_{ij}a_{jl}c_l$; exact: $\tfrac1{24}$): $ \text{(O8)}\qquad b_4a_{43}a_{32}c_2=\frac1{24}. $ **What has been established.** Exactly eight independent conditions (O1)–(O8) are necessary and sufficient for the four-stage explicit method (with consistency relations) to reproduce the exact Taylor expansion through $h^4$ for every $C^4$ function $f$. Each condition has been traced to the unique elementary-differential group that generates it, and the potential double-counting from overlapping monomials has been shown to be vacuous. --- # 8. The order conditions in scalar and Butcher form ## 8.1 Expanded scalar form (four stages, $c_1=0$) $ \begin{aligned} &\text{(O1)} & b_1+b_2+b_3+b_4&=1\\ &\text{(O2)} & b_2c_2+b_3c_3+b_4c_4&=\tfrac12\\ &\text{(O3)} & b_2c_2^2+b_3c_3^2+b_4c_4^2&=\tfrac13\\ &\text{(O4)} & b_3a_{32}c_2+b_4a_{42}c_2+b_4a_{43}c_3&=\tfrac16\\ &\text{(O5)} & b_2c_2^3+b_3c_3^3+b_4c_4^3&=\tfrac14\\ &\text{(O6)} & b_3c_3a_{32}c_2+b_4c_4a_{42}c_2+b_4c_4a_{43}c_3&=\tfrac18\\ &\text{(O7)} & b_3a_{32}c_2^2+b_4a_{42}c_2^2+b_4a_{43}c_3^2&=\tfrac1{12}\\ &\text{(O8)} & b_4a_{43}a_{32}c_2&=\tfrac1{24} \end{aligned} $ ## 8.2 Compact Butcher form Let $A=(a_{ij})\in\mathbb{R}^{4\times4}$ (strictly lower triangular), $b=(b_1,\dots,b_4)^{\mathsf T}$, $c=(0,c_2,c_3,c_4)^{\mathsf T}$, $\mathbf e=(1,1,1,1)^{\mathsf T}$, and let $c^2,c^3$ denote componentwise powers and $c\circ v$ the componentwise (Hadamard) product. The consistency relations read $A\mathbf e=c$, and the eight conditions become $ \begin{array}{llll} \text{(O1)}\;\; b^{\mathsf T}\mathbf e=1, &\text{(O2)}\;\; b^{\mathsf T}c=\dfrac12, &\text{(O3)}\;\; b^{\mathsf T}c^2=\dfrac13, &\text{(O4)}\;\; b^{\mathsf T}Ac=\dfrac16,\\[8pt] \text{(O5)}\;\; b^{\mathsf T}c^3=\dfrac14, &\text{(O6)}\;\; b^{\mathsf T}(c\circ Ac)=\dfrac18, &\text{(O7)}\;\; b^{\mathsf T}Ac^2=\dfrac1{12}, &\text{(O8)}\;\; b^{\mathsf T}A^2c=\dfrac1{24}. \end{array} $ **What has been established.** The complete order-4 condition set in both notations. The Butcher forms are transliterations of (O1)–(O8), no more and no less: e.g. $b^{\mathsf T}A^2c=\sum_{i,j,l}b_ia_{ij}a_{jl}c_l$ reduces, for a strictly lower-triangular $A$ with $c_1=0$, to the single product $b_4a_{43}a_{32}c_2$ appearing in (O8). --- # 9. The structure of the system: polynomial and nonlinear The unknowns are the thirteen numbers $ c_2,c_3,c_4,\;a_{21},a_{31},a_{32},a_{41},a_{42},a_{43},\;b_1,b_2,b_3,b_4, $ subject to the three consistency relations (which eliminate $a_{21}$ and one entry from each later row, leaving ten free parameters) and the eight conditions (O1)–(O8). Each condition is a *polynomial* equation in the tableau entries — inevitably so, because the stage expansions are compositions of polynomials in those entries — but the system is not linear, and it is worth being precise about where the nonlinearity lives. **Linear in $b$ once $(A,c)$ is fixed.** Every condition (O1)–(O8) has the form $\sum_i b_i\,\varphi_i(A,c)=\text{const}$: each is *linear in $b$* with coefficients $\varphi_i$ that are polynomials in $A$ and $c$. Explicitly, the eight functionals are $ \varphi\in\Bigl\{1,\;c_i,\;c_i^2,\;\textstyle\sum_ja_{ij}c_j,\;c_i^3,\;c_i\sum_ja_{ij}c_j,\;\sum_ja_{ij}c_j^2,\;\sum_{j,l}a_{ij}a_{jl}c_l\Bigr\}. $ So for *fixed* $(A,c)$, finding $b$ means solving eight linear equations in four unknowns — generically unsolvable, which is the correct intuition for why the nodes and the coupling matrix cannot be chosen arbitrarily: $(A,c)$ must be tuned so that the eight linear functionals become compatible (they must, in effect, collapse onto a four-dimensional solvable system). **Nonlinear across parameter blocks.** The nonlinearity is of two kinds: - *Products of $b$ with $(A,c)$:* every condition beyond (O1) multiplies $b_i$ against powers of $c$ and entries of $A$; e.g. (O2) is bilinear in $(b,c)$, and (O6) is trilinear (degree one in each of $b$, $c$, and the pair $(A,c)$ jointly through $a_{ij}c_j$). - *Products internal to $(A,c)$:* condition (O8), $b_4a_{43}a_{32}c_2=\tfrac1{24}$, is genuinely cubic in the coupling entries and quartic overall; it cannot be linearized by fixing any single block. It is the algebraic fingerprint of double nesting (stage 4 reading stage 3 reading stage 2), just as (O4) with its single product $a_{ij}c_j$ fingerprints single nesting. Note also the internal hierarchy: (O1), (O2), (O3), (O5) involve only $(b,c)$ — they say that the pairs $(b_i,c_i)$ form a *quadrature rule* on $[0,1]$ exact for polynomials $1,x,x^2,x^3$: $ \sum_i b_ic_i^{q}=\frac1{q+1},\qquad q=0,1,2,3, $ i.e. exact for all cubics. The remaining four conditions (O4), (O6), (O7), (O8) constrain how $A$ threads the stages together. This split organizes the solution in §10. **What has been established.** The system (O1)–(O8) is polynomial of total degree up to four; it is linear in $b$ for fixed $(A,c)$; the quadrature subsystem involves $(b,c)$ alone; and the coupling conditions — above all the doubly nested (O8) — are irreducibly nonlinear in the tableau as a whole. --- # 10. Deriving the classical coefficients We now *construct* a solution by a short sequence of explicitly flagged choices, solving for everything else. Choices are marked **[choice]**; consequences forced by the order conditions are marked **[forced]**. ## 10.1 [choice] The quadrature skeleton: $c_2=c_3=\tfrac12$, $c_4=1$ Motivation. By §9 the weights and nodes must form a quadrature rule exact on cubics. The lowest-cardinality classical rule with that exactness on $[0,1]$ using the endpoints is Simpson's rule, with nodes $0,\tfrac12,1$ and weights $\tfrac16,\tfrac46,\tfrac16$. A four-stage method has four nodes with $c_1=0$; the most economical way to reach the Simpson pattern is to *double the midpoint node*, $ c=(0,\tfrac12,\tfrac12,1)^{\mathsf T}, $ letting the two midpoint stages share Simpson's middle weight. Nothing yet dictates how that weight splits between $b_2$ and $b_3$. (This choice is where the "classical" in classical RK4 enters; other node vectors, e.g. $c=(0,\tfrac13,\tfrac23,1)$, lead to the equally valid 3/8-rule method.) The consistency relation for row two is then **[forced]**: $ a_{21}=c_2=\tfrac12. $ ## 10.2 [forced] The weights, from the quadrature subsystem With the chosen nodes, (O1), (O2), (O3), (O5) read $ \begin{aligned} b_1+b_2+b_3+b_4&=1,\\ \tfrac12(b_2+b_3)+b_4&=\tfrac12,\\ \tfrac14(b_2+b_3)+b_4&=\tfrac13,\\ \tfrac18(b_2+b_3)+b_4&=\tfrac14. \end{aligned} $ Only the combination $s:=b_2+b_3$ enters. Subtract the third equation from the second: $ \Bigl(\tfrac12-\tfrac14\Bigr)s=\tfrac12-\tfrac13 \;\Longrightarrow\; \tfrac14 s=\tfrac16 \;\Longrightarrow\; s=\tfrac23. $ Then from the second equation $b_4=\tfrac12-\tfrac12\cdot\tfrac23=\tfrac12-\tfrac13=\tfrac16$, and from the first $b_1=1-\tfrac23-\tfrac16=\tfrac16$. The fourth equation is now a consistency check, not a constraint: $ \tfrac18\cdot\tfrac23+\tfrac16=\tfrac1{12}+\tfrac2{12}=\tfrac3{12}=\tfrac14.\;\checkmark $ (It had to check: Simpson's rule is exact on cubics, and node-doubling preserves exactness whatever the split of the middle weight.) So $ b_1=\tfrac16,\qquad b_2+b_3=\tfrac23,\qquad b_4=\tfrac16, $ with the split of $\tfrac23$ still free. ## 10.3 [forced] Consequences of the coupling conditions Insert $c_2=c_3=\tfrac12$, $c_4=1$, $b_4=\tfrac16$ into (O4), (O6), (O7), (O8). Write $u:=a_{32}$ and $w:=a_{42}+a_{43}$ for brevity. $ \text{(O4)}:\quad b_3u\cdot\tfrac12+\tfrac16\Bigl(a_{42}\cdot\tfrac12+a_{43}\cdot\tfrac12\Bigr)=\tfrac16 \;\Longrightarrow\; \tfrac12 b_3u+\tfrac1{12}w=\tfrac16. $ $ \text{(O6)}:\quad b_3\cdot\tfrac12\cdot u\cdot\tfrac12+\tfrac16\cdot1\cdot\Bigl(a_{42}\cdot\tfrac12+a_{43}\cdot\tfrac12\Bigr)=\tfrac18 \;\Longrightarrow\; \tfrac14 b_3u+\tfrac1{12}w=\tfrac18. $ $ \text{(O7)}:\quad b_3u\cdot\tfrac14+\tfrac16\Bigl(a_{42}\cdot\tfrac14+a_{43}\cdot\tfrac14\Bigr)=\tfrac1{12} \;\Longrightarrow\; \tfrac14 b_3u+\tfrac1{24}w=\tfrac1{12}. $ $ \text{(O8)}:\quad \tfrac16\,a_{43}\,u\cdot\tfrac12=\tfrac1{24} \;\Longrightarrow\; a_{43}u=\tfrac12. $ Subtract (O7)′ from (O6)′: $ \Bigl(\tfrac1{12}-\tfrac1{24}\Bigr)w=\tfrac18-\tfrac1{12} \;\Longrightarrow\; \tfrac1{24}w=\tfrac1{24} \;\Longrightarrow\; w=a_{42}+a_{43}=1. $ Since $c_4=a_{41}+a_{42}+a_{43}=1$, this **forces** $ a_{41}=0 $ — the fourth stage cannot lean on $k_1$ at all; this is not a stylistic choice but a theorem, given the node vector. Substituting $w=1$ back into (O6)′: $ \tfrac14b_3u=\tfrac18-\tfrac1{12}=\tfrac1{24} \;\Longrightarrow\; b_3u=\tfrac16, $ and then (O4)′ is automatically satisfied: $\tfrac12\cdot\tfrac16+\tfrac1{12}=\tfrac1{12}+\tfrac1{12}=\tfrac16.\;\checkmark$ At this point the entire residual freedom (for the chosen nodes) is the single parameter $u=a_{32}\neq0$: $ b_3=\frac{1}{6u},\qquad b_2=\frac23-\frac1{6u},\qquad a_{43}=\frac1{2u},\qquad a_{42}=1-\frac1{2u},\qquad a_{31}=\frac12-u, $ a one-parameter family of fourth-order methods all sharing $c=(0,\tfrac12,\tfrac12,1)$. ## 10.4 [choice] Closing the family: $a_{32}=\tfrac12$ The classical construction takes $ u=a_{32}=\tfrac12. $ Motivation, in decreasing order of force: (i) it makes $a_{31}=\tfrac12-u=0$, so stage 3 is a *pure midpoint correction* — a half-step taken with the freshest available slope $k_2$ only, mirroring how stage 2 took a half-step with $k_1$ only; every stage then uses exactly one previous slope, the most recent one. (ii) It symmetrizes the weights: $b_3=\tfrac1{6u}=\tfrac13$ and $b_2=\tfrac23-\tfrac13=\tfrac13$, splitting Simpson's middle weight evenly across the two midpoint stages. (iii) It yields $a_{43}=\tfrac1{2u}=1$ and $a_{42}=0$, so stage 4 is a *full* step taken with the freshest slope $k_3$ only. The resulting evaluation pattern — Euler half-step, midpoint half-step, midpoint full-step, Simpson combination — is the memorable "slope at the start, two slopes at the middle, slope at the end." **[forced]** With $u=\tfrac12$: $ b_2=b_3=\tfrac13,\qquad a_{31}=0,\quad a_{32}=\tfrac12,\qquad a_{41}=0,\quad a_{42}=0,\quad a_{43}=1. $ ## 10.5 Inventory of choices versus consequences - **Chosen:** $c_2=c_3=\tfrac12$, $c_4=1$ (Simpson skeleton with doubled midpoint); $a_{32}=\tfrac12$ (fresh-slope/symmetry normalization). - **Forced by consistency:** $a_{21}=\tfrac12$; $a_{31}=c_3-a_{32}$; $a_{41}=c_4-a_{42}-a_{43}$. - **Forced by the order conditions:** $b_1=b_4=\tfrac16$; $b_2+b_3=\tfrac23$; $a_{41}=0$; $b_3a_{32}=\tfrac16$; $a_{43}a_{32}=\tfrac12$ — whence, after the single choice $a_{32}=\tfrac12$: $b_2=b_3=\tfrac13$, $a_{31}=0$, $a_{42}=0$, $a_{43}=1$. **What has been established.** The classical tableau is not verified but *derived*: two flagged choices (the node vector and the value of $a_{32}$) reduce the nonlinear system to forced linear consequences, and the residual one-parameter family containing the classical method has been exhibited explicitly. --- # 11. The classical tableau and direct verification of every condition $ \begin{array}{c|cccc} 0 & 0 & 0 & 0 & 0\\[2pt] \tfrac12 & \tfrac12 & 0 & 0 & 0\\[2pt] \tfrac12 & 0 & \tfrac12 & 0 & 0\\[2pt] 1 & 0 & 0 & 1 & 0\\[2pt] \hline & \tfrac16 & \tfrac13 & \tfrac13 & \tfrac16 \end{array} $ i.e. $b=\bigl(\tfrac16,\tfrac13,\tfrac13,\tfrac16\bigr)^{\mathsf T}$, $c=\bigl(0,\tfrac12,\tfrac12,1\bigr)^{\mathsf T}$, and the stages $ k_1=f(t,y),\quad k_2=f\Bigl(t+\tfrac h2,\,y+\tfrac h2k_1\Bigr),\quad k_3=f\Bigl(t+\tfrac h2,\,y+\tfrac h2k_2\Bigr),\quad k_4=f\bigl(t+h,\,y+hk_3\bigr), $ $ y_{n+1}=y+\frac h6\bigl(k_1+2k_2+2k_3+k_4\bigr). $ Verification, one condition at a time: $ \text{(O1)}:\quad \tfrac16+\tfrac13+\tfrac13+\tfrac16=\tfrac{1+2+2+1}{6}=1.\;\checkmark $ $ \text{(O2)}:\quad \tfrac13\cdot\tfrac12+\tfrac13\cdot\tfrac12+\tfrac16\cdot1 =\tfrac16+\tfrac16+\tfrac16=\tfrac12.\;\checkmark $ $ \text{(O3)}:\quad \tfrac13\cdot\tfrac14+\tfrac13\cdot\tfrac14+\tfrac16\cdot1 =\tfrac1{12}+\tfrac1{12}+\tfrac2{12}=\tfrac{4}{12}=\tfrac13.\;\checkmark $ $ \text{(O4)}:\quad b_3a_{32}c_2+b_4a_{43}c_3 =\tfrac13\cdot\tfrac12\cdot\tfrac12+\tfrac16\cdot1\cdot\tfrac12 =\tfrac1{12}+\tfrac1{12}=\tfrac16.\;\checkmark $ (The $b_4a_{42}c_2$ term vanishes since $a_{42}=0$.) $ \text{(O5)}:\quad \tfrac13\cdot\tfrac18+\tfrac13\cdot\tfrac18+\tfrac16\cdot1 =\tfrac1{24}+\tfrac1{24}+\tfrac4{24}=\tfrac{6}{24}=\tfrac14.\;\checkmark $ $ \text{(O6)}:\quad b_3c_3a_{32}c_2+b_4c_4a_{43}c_3 =\tfrac13\cdot\tfrac12\cdot\tfrac12\cdot\tfrac12+\tfrac16\cdot1\cdot1\cdot\tfrac12 =\tfrac1{24}+\tfrac1{12}=\tfrac{1+2}{24}=\tfrac18.\;\checkmark $ $ \text{(O7)}:\quad b_3a_{32}c_2^2+b_4a_{43}c_3^2 =\tfrac13\cdot\tfrac12\cdot\tfrac14+\tfrac16\cdot1\cdot\tfrac14 =\tfrac1{24}+\tfrac1{24}=\tfrac1{12}.\;\checkmark $ $ \text{(O8)}:\quad b_4a_{43}a_{32}c_2=\tfrac16\cdot1\cdot\tfrac12\cdot\tfrac12=\tfrac1{24}.\;\checkmark $ **What has been established.** The classical tableau satisfies all eight conditions (O1)–(O8); by §7 its one-step expansion therefore agrees with the exact solution's Taylor expansion through $h^4$ for every $f\in C^4$. --- # 12. Local truncation error $O(h^5)$ and global error $O(h^4)$ ## 12.1 The increment function and the local truncation error Write the method as a one-step map with *increment function* $ \Phi(t,y;h):=\sum_{i=1}^4 b_i\,k_i(t,y;h), $ where the $k_i(t,y;h)$ are the stage values computed from base point $(t,y)$; note $\Phi$ is defined for any $(t,y)$, not only points on the exact solution. The *local truncation error* at step $n$ is the defect of the exact solution in the scheme: $ \tau_{n+1}:=y(t_n+h)-\Bigl[y(t_n)+h\,\Phi\bigl(t_n,y(t_n);h\bigr)\Bigr]. $ **Claim: $\tau_{n+1}=O(h^5)$, uniformly for $t_n$ in a compact interval.** Both quantities inside the bracket were expanded about $(t_n,y(t_n))$ in §§1 and 6. By construction the coefficients of $h^0,h^1,h^2,h^3,h^4$ — that is, of every one of the 19 monomials of §2 in every power — agree, because (O1)–(O8) hold (§11). Therefore $\tau_{n+1}$ equals the difference of the two remainder terms. Each remainder is explicitly controlled: the exact-solution remainder is $\tfrac{h^5}{120}y^{(5)}(\xi)$, and $y^{(5)}$ is a polynomial in the partials of $f$ through total order four evaluated along the solution — continuous, hence bounded on a compact tube $ K=\bigl\{(t,z): t_0\le t\le T,\ |z-y(t)|\le 1\bigr\} $ around the solution graph. The method-side remainder is a finite sum of terms of the form (tableau polynomial) $\times$ (fourth-order Taylor remainders $R_4$ of §3.1, and products of lower remainders with bounded factors), each bounded by a constant times $h^5$ after the overall factor $h$ from the update, with constants again controlled by $\sup_K$ of the partials of $f$ through order four. Hence there exist $C>0$ and $h_0>0$, depending only on those suprema and the tableau, with $ |\tau_{n+1}|\le C\,h^5 \qquad\text{for all }0<h\le h_0,\ t_n\in[t_0,T-h]. $ This is the promised fifth-order local error; the method is *consistent of order four*. ## 12.2 Stability of the increment function Assume, as standard, that $f$ is Lipschitz in $y$ on $K$, uniformly in $t$: $ |f(t,y)-f(t,z)|\le L|y-z| \qquad\text{for }(t,y),(t,z)\in K. $ Then each stage is Lipschitz in the base value $y$, by induction along the stages. Indeed $|k_1(t,y)-k_1(t,z)|\le L|y-z|=:L_1|y-z|$, and if $|k_j(t,y)-k_j(t,z)|\le L_j|y-z|$ for $j<i$, then $ |k_i(t,y)-k_i(t,z)| =\Bigl|f\Bigl(t+c_ih,\,y+h\sum_{j<i}a_{ij}k_j(t,y)\Bigr)-f\Bigl(t+c_ih,\,z+h\sum_{j<i}a_{ij}k_j(t,z)\Bigr)\Bigr| \le L\Bigl(1+h\sum_{j<i}|a_{ij}|L_j\Bigr)|y-z|, $ so $L_i:=L\bigl(1+h\sum_{j<i}|a_{ij}|L_j\bigr)$ works. For $0<h\le h_0$ the $L_i$ are bounded (they are polynomials in $hL$ with tableau coefficients; for the classical tableau, $L_2=L(1+\tfrac{hL}{2})$, $L_3=L(1+\tfrac{h}{2}L_2)$, $L_4=L(1+hL_3)$), and consequently $ |\Phi(t,y;h)-\Phi(t,z;h)|\le\Lambda|y-z|, \qquad \Lambda:=\sum_{i=1}^4|b_i|\,L_i\Big|_{h=h_0} $ (the $L_i$ are increasing in $h$, so evaluating at $h_0$ gives a uniform constant). This Lipschitz property of $\Phi$ is the *stability* half of the convergence argument. ## 12.3 From local to global: the error recursion Let $y_0=y(t_0)$ (exact start), $t_n=t_0+nh$, $N=\lfloor (T-t_0)/h\rfloor$, and $e_n:=|y(t_n)-y_n|$. Insert-and-subtract using the definition of $\tau_{n+1}$: $ \begin{aligned} y(t_{n+1})-y_{n+1} &=\underbrace{y(t_{n+1})-y(t_n)-h\Phi\bigl(t_n,y(t_n);h\bigr)}_{=\;\tau_{n+1}} +\underbrace{y(t_n)-y_n}_{\text{propagated error}} +h\Bigl[\Phi\bigl(t_n,y(t_n);h\bigr)-\Phi\bigl(t_n,y_n;h\bigr)\Bigr]. \end{aligned} $ Taking absolute values and applying §12.1 and §12.2 (valid as long as the numerical trajectory stays in $K$, which holds for $h$ small enough by a short bootstrap argument: the bound below itself keeps $e_n\le1$ for $h\le h_1\le h_0$): $ e_{n+1}\le(1+h\Lambda)\,e_n+C h^5. $ Unroll the recursion with $e_0=0$. One application gives $e_1\le Ch^5$; two give $e_2\le(1+h\Lambda)Ch^5+Ch^5$; inductively, $ e_n\le Ch^5\sum_{j=0}^{n-1}(1+h\Lambda)^j =Ch^5\,\frac{(1+h\Lambda)^n-1}{(1+h\Lambda)-1} =\frac{Ch^4}{\Lambda}\Bigl[(1+h\Lambda)^n-1\Bigr], $ by the geometric-series formula. Finally, from $1+x\le e^x$ applied with $x=h\Lambda$, $ (1+h\Lambda)^n\le e^{nh\Lambda}\le e^{\Lambda(T-t_0)} \qquad(nh\le T-t_0), $ so that $ \boxed{\; \max_{0\le n\le N}\;|y(t_n)-y_n| \;\le\; \frac{C\bigl(e^{\Lambda(T-t_0)}-1\bigr)}{\Lambda}\;h^4 \;=\;O(h^4).\;} $ The mechanism deserves emphasis: the method commits $N\approx(T-t_0)/h$ local errors of size $O(h^5)$; stability lets them accumulate at worst geometrically with ratio $1+h\Lambda$ per step, and the sum of $O(1/h)$ terms of size $O(h^5)$, amplified by the bounded factor $e^{\Lambda(T-t_0)}$, is $O(h^4)$. One power of $h$ is always lost passing from local to global; a fifth-order local error is exactly what a fourth-order method means. **What has been established.** For $f\in C^4$ and Lipschitz in $y$ near the solution: the classical RK4 one-step error is $O(h^5)$ with an explicit constant-generating argument, the increment function is Lipschitz with a computable constant, and the discrete Grönwall/telescoping argument converts these into the global bound $O(h^4)$ on any fixed compact interval. --- # Postscript: the rooted-tree ledger (after the fact) Only now, with the calculation complete, is it worth recording the combinatorial pattern it instantiates. Each of the eight groups matched in §7 corresponds to a rooted tree with at most four vertices — for the scalar nonautonomous problem, each vertex represents an $f$-evaluation, and a vertex with $m$ children contributes a mixed partial of total order $m$ distributed between $t$ and $y$ (which is why each "tree" showed up here as a *bundle* of monomials: the group $G$, for instance, is the bushy three-vertex tree, its three monomials $f_{tt}, ff_{ty}, f^2f_{yy}$ being the binomial split of a second directional derivative). The dictionary for our eight conditions: | Group | Tree (bracket notation) | Condition | Value $=1/\gamma$ | |---|---|---|---| | $f$ | $\bullet$ | $b^{\mathsf T}\mathbf e$ | $1$ | | $F$ | $[\bullet]$ | $b^{\mathsf T}c$ | $1/2$ | | $G$ | $[\bullet,\bullet]$ | $b^{\mathsf T}c^2$ | $1/3$ | | $f_yF$ | $[[\bullet]]$ | $b^{\mathsf T}Ac$ | $1/6$ | | $H$ | $[\bullet,\bullet,\bullet]$ | $b^{\mathsf T}c^3$ | $1/4$ | | $FP$ | $[\bullet,[\bullet]]$ | $b^{\mathsf T}(c\circ Ac)$ | $1/8$ | | $f_yG$ | $[[\bullet,\bullet]]$ | $b^{\mathsf T}Ac^2$ | $1/12$ | | $f_y^2F$ | $[[[\bullet]]]$ | $b^{\mathsf T}A^2c$ | $1/24$ | The right-hand sides $1,\tfrac12,\tfrac13,\tfrac16,\tfrac14,\tfrac18,\tfrac1{12},\tfrac1{24}$ are the reciprocals of the tree *densities* $\gamma(\tau)$, and the exact-side coefficients $\tfrac1{q!}\times(\text{multiplicity})$ that we computed by hand in §1.4 are the statement $\alpha(\tau)=$ (number of monotone labelings)$/\,$symmetry, specialized to one dimension. There are $1,1,2,4$ rooted trees of orders $1,2,3,4$ — hence eight conditions — and the count grows to $9$ additional trees at order five, which is the beginning of the explanation for why fifth order requires six stages rather than five. None of this was needed above; it is the bookkeeping that makes order six and beyond humane. --- # Closing summary Starting only from the explicit four-stage ansatz and the smoothness of $f$: (1) the exact flow expands as $y+hf+\tfrac{h^2}{2}F+\tfrac{h^3}{6}(G+f_yF)+\tfrac{h^4}{24}(H+f_yG+3FP+f_y^2F)+O(h^5)$; (2) each stage, under the internally derived consistency relations $c_i=\sum_ja_{ij}$, expands in the same eight elementary-differential groups with tableau-polynomial coefficients; (3) matching forces exactly the eight conditions (O1)–(O8); (4) choosing the doubled-midpoint Simpson skeleton $c=(0,\tfrac12,\tfrac12,1)$ and the fresh-slope normalization $a_{32}=\tfrac12$ forces every remaining entry, yielding the classical tableau with weights $\tfrac16,\tfrac13,\tfrac13,\tfrac16$; (5) the resulting one-step error is $O(h^5)$ and the global error $O(h^4)$ by stability plus telescoping. The classical RK4 method is, in this precise sense, Simpson's rule taught to read a slope field.