Integration, accumulation, and simple differential equations

A derivative describes a local rate. An integral accumulates a rate across an interval, with units, signs and limiting hypotheses that a finite sum alone does not provide. A differential equation then specifies how a state changes through time. Its numerical simulation is a separate approximation, and a stable mathematical solution can be simulated badly by an unsuitable step.

Retrieval check: recall Module16 product/chain rules, mean-value control and Taylor error, plus Module15 completeness and continuity. The labs need only standard-library Python. Density examples introduce the integral normalisation contract; probability models and expectations are developed systematically in Modules21–23.

Partition a closed interval[a,b], a<b, into a=x₀<x₁<…<x_n=b. Each subinterval width is Δx_i=x_i−x_{i−1}>0, and choose a tag ξ_i inside it. A Riemann sum is Σf(ξ_i)Δx_i. The mesh is the largest width, not merely the number of sample points. The definite integral is a common finite limit approached by all such tagged sums as the mesh tends to zero. A specific uniform grid can illustrate this process, but agreement along one chosen sampling pattern alone is weaker than the integrability definition.

For a bounded function, let lower and upper sums use each subinterval’s infimum and supremum instead of a tag. Every tagged sum lies between them. Refining a partition increases its lower sum and decreases its upper sum. Lower sums from one partition do not exceed upper sums from another, because a common refinement allows comparison. Completeness therefore supplies the supremum of all lower sums and infimum of all upper sums. If they agree, the function is Riemann integrable and that shared value is its integral. This is a useful equivalent construction of the limit rather than a numerical integration recipe by itself.

A continuous function on a closed bounded interval is uniformly continuous: one input tolerance controls output differences for every pair throughout the interval. To see why pointwise continuity suffices here, suppose uniform continuity failed. Then pairs at distances tending to zero could have output distances at least some fixed ε. Compactness supplies a convergent subsequence of the first inputs; the second inputs share its limit, and continuity makes their output difference vanish, a contradiction. The finite interval hypotheses are doing the work. Pointwise continuity on an open or unbounded domain need not provide the same uniform guarantee.

Uniform continuity proves integrability. Choose a mesh small enough that all values within any subinterval differ by less than ε. Its upper-minus-lower sum is then at most εΣΔx_i=ε(b−a). Such gaps can be made arbitrarily small, so the two completeness bounds agree. It also shows arbitrary sufficiently fine tagged sums approach the same number. An integral exists for other functions too, including bounded piecewise-continuous functions with finitely many jumps, but not every bounded function is Riemann integrable. A function equal to one on rational points and zero on irrational points has every lower sum zero and every upper sum b−a.

The integral is signed accumulation. Positive portions add, negative portions subtract. Total geometric area uses ∫|f| instead. For f(x)=x on[−1,1], signed integral is zero while total area is one. Define ∫_a^a f=0 and reversed limits by ∫_b^a f=−∫_a^b f. Linearity, interval additivity and order preservation follow first for finite sums and then by limits. In particular, if m≤f≤M on[a,b], then m(b−a)≤∫f≤M(b−a), and |∫f|≤∫|f|. Output units multiply input units: metres per second integrated over seconds gives metres.

For continuous f on[a,b], define accumulation A(x)=∫_a^x f(t)dt. At an interior x, A(x+h)−A(x)=∫_x^{x+h}f(t)dt. Subtract h f(x), divide by h and use continuity: all nearby f(t) differ from f(x) by at most any requested ε, so the absolute quotient error is at most ε. The argument works for either sign of h using reversed orientation. Therefore A′(x)=f(x). This is the first fundamental-theorem connection: a continuous rate is the derivative of its accumulated quantity. The integration variable t is a bound dummy variable; it is different from the moving endpoint x.

An antiderivative F satisfies F′=f on an interval. Since A is another such function, (F−A)′=0; the mean value theorem makes the difference constant. Evaluating at endpoints gives the second connection, ∫_a^b f=F(b)−F(a). This proves the familiar evaluation rule under the stated continuity and derivative conditions. A definite integral is one number for fixed endpoints. An indefinite integral is notation for a family F+C, where C is arbitrary; it is not a fixed accumulated number until a boundary or initial condition selects the constant.

Basic antiderivatives follow by reversing verified derivative formulas. On appropriate intervals, x^p integrates to x^{p+1}/(p+1) when p≠−1; use positive x for general real powers, and broader domains for polynomial integer cases. The excluded p=−1 gives ln|x| on each interval not crossing zero. exp integrates to exp, cosine to sine, and sine to minus cosine. Linear combinations follow by linearity. A table does not eliminate domain checks, especially when an antiderivative expression exists only on one connected component.

Constants can differ across disconnected domains. If f(x)=1/x on R excluding zero, every antiderivative has form ln|x|+C₋ on the negative component and ln|x|+C₊ on the positive one. A derivative does not constrain the difference of those constants across the missing point. Similarly, writing F(b)−F(a) across a singularity is not automatically a proper integral; the fundamental theorem requires an appropriate whole interval. The improper limits in Section4 analyse those cases explicitly.

Moving endpoints add another dependency layer. If f is continuous and u,v are differentiable with their ranges inside the integration interval, then B(x)=∫_{u(x)}^{v(x)}f(t)dt=A(v(x))−A(u(x)), so B′(x)=f(v(x))v′(x)−f(u(x))u′(x). Both endpoints contribute, with opposite signs. For ∫_0^{x²}exp(t)dt, the derivative is2x exp(x²), not just exp(x²). This is the fundamental theorem plus the chain rule, not cancellation of formal symbols without hypotheses.

The first theorem above assumes continuous f. For discontinuous but integrable f, accumulation can still be continuous, while its derivative at a jump may fail or need separate analysis. At a point of continuity the local average argument remains valid. Do not assert A′ equals an arbitrarily assigned point value at every discontinuity. Changing a bounded integrand at finitely many points leaves its integral unchanged, because those points can be enclosed in intervals of arbitrarily small total width. The accumulation process therefore need not encode those isolated assignments in a derivative everywhere.

Substitution follows from the chain rule and fundamental theorem. If g is continuously differentiable on[a,b] and f is continuous on an interval containing its range, choose F′=f. Then [F(g(x))]′=f(g(x))g′(x), so ∫a^b f(g(x))g′(x)dx=F(g(b))−F(g(a))=∫{g(a)}^{g(b)}f(u)du. For this form, g need not be monotone; the endpoint expression automatically records orientation and any cancellation along the path. An interpretation that inverts g or removes its derivative factor may need additional injectivity or branch assumptions.

Not every resemblance permits a substitution. The expression exp(x²) alone lacks the factor2x, so the preceding elementary cancellation does not evaluate its integral. A substitution can still change variables, but it may introduce a new factor or non-elementary integral rather than simplify it. Likewise, taking u=x² on a domain crossing zero and trying to write dx=du/(2sqrt(u)) without separating branches loses the negative-x branch. The chain-rule identity remains sound, while an inverse-variable derivation has a different set of domain requirements.

Integration by parts follows from (uv)′=u′v+uv′. Integrate this identity over[a,b] and rearrange: ∫_a^b u(x)v′(x)dx=[u(x)v(x)]_a^b−∫_a^b u′(x)v(x)dx, for continuously differentiable u,v. One factor is differentiated and the other must be integrated, so choices matter. The boundary term is part of the formula; omitting it changes the definite result even when the remaining integrand looks simpler. Indefinite notation packages the same identity plus an arbitrary constant.

For ∫_0^1 x exp(x)dx, choose u=x and v′=exp(x), hence u′=1,v=exp(x). The boundary term is e, and the remaining integral is e−1, giving result one. Reversing the choice would require an antiderivative of x and leave a different residual integral. The useful choice is usually the one reducing complexity, but success is mathematical rather than guaranteed by a mnemonic. Repeated integration by parts can handle polynomial times exponential terms by lowering polynomial degree.

For ∫_1^e ln x dx, choose u=ln x,v′=1, so v=x. The boundary term[x ln x]_1^e=e and ∫_1^e x(1/x)dx=e−1, again giving one. The logarithm domain is positive throughout. For an improper integral, apply integration by parts on truncated proper intervals first and then check each boundary and residual limit. A formal identity with a divergent boundary expression is not a finite evaluation.

These transformations can be checked by differentiation for indefinite claims and by explicit boundaries for definite ones. Differentiating a proposed antiderivative is often simpler than finding it, but the check must happen on its claimed domain. Multiple antiderivative forms may differ by a constant on a connected interval. If their difference is not constant or a domain piece has been lost, the formulas do not describe the same family there. In applications, units also provide a useful check: substitution changes the input coordinate, while the derivative factor accounts for its rate of scaling.

A proper integral of a continuous function uses a finite closed interval without an unbounded integrand. Improper integrals extend that definition through limits. At a singular left endpoint, ∫a^b f means lim{ε↓0}∫_{a+ε}^b f, if the finite limit exists. At infinity, ∫a^∞ f means lim{R→∞}∫_a^R f. A singularity inside an interval requires splitting into left and right improper integrals; each must converge separately. Two divergent sides cannot be cancelled and called an ordinary convergent improper integral.

For example ∫_ε^1 1/x dx=−ln ε grows without bound as ε↓0, so ∫_0^1 1/x diverges. Every finite cutoff yields a finite number, but the limiting claim fails. In contrast, ∫_ε^1 x^{−1/2}dx=2(1−sqrt(ε))→2, so an unbounded integrand can have a finite improper integral. Bounded height is not necessary for this extension, and a finite cutoff calculation is not sufficient for convergence. The nature of the singularity matters.

More generally, for p≠1, ∫_ε^1 x^{−p}dx=(1−ε^{1−p})/(1−p). This has a finite limit exactly when p<1; p=1 is the divergent logarithm case. On[1,∞), ∫_1^R x^{−p}dx=(R^{1−p}−1)/(1−p), which converges exactly when p>1. The same exponent family has opposite endpoint and infinity criteria. Both statements require the corresponding limiting calculation, rather than an informal assertion that a small function or an unbounded function must behave one way.

The signed example1/x on[−1,1] has separately divergent left and right integrals. Symmetric truncation ∫_{−1}^{−ε}1/x+∫_ε^1 1/x=0 for every ε gives a Cauchy principal value of zero, a separately defined symmetric-limit quantity. It does not make either side converge or turn the ordinary improper integral into zero. Keeping different notions of limit distinct avoids an apparent cancellation proof of a nonexistent normalisation constant.

A probability-density candidate must be nonnegative and have total integral one over its whole support. For an unnormalised nonnegative q with finite positive integral Z, define p=q/Z. If Z=0 or infinite, that division does not produce a valid unit-mass density. Sign-changing q can have integral one while still failing nonnegativity. A grid’s sampled nonnegative values do not certify that condition everywhere without structural evidence. In the linear lab, c≥0 proves it for all x∈[0,1]; the proof is not replaced by the grid.

For λ>0, q(x)=exp(−λx) on x≥0 has integral1/λ, so p(x)=λexp(−λx). Its mass on[0,t] is1−exp(−λt). Integration by parts on[0,R] gives the first moment before taking R→∞; the boundary R exp(−λR) vanishes, yielding mean1/λ. The vanishing boundary can be justified using exp(λR)≥(λR)²/2 from its positive series, so R exp(−λR)≤2/(λ²R)→0. This is a complete limit check, not just writing a zero at infinity.

Probability models will clarify event definitions and expectations later. Here the mathematical distinction already matters: probability assigned to an interval is an integral, while a single point has zero mass for these continuous density models. Changing a density at one isolated point does not change its integrals. Density units are reciprocal input units, so its numerical height can change under a unit conversion without changing event probabilities. Keep the support, scale and integration measure explicit.

Quadrature approximates a definite integral with a weighted finite sum. For n equal subintervals, h=(b−a)/n and x_i=a+ih. The composite trapezoidal rule uses h[f(a)/2+Σ_{i=1}^{n−1}f(x_i)+f(b)/2]. It integrates the piecewise linear interpolant through sampled values, giving endpoint half-weights because adjacent trapezoids share interior values. It is exact for linear functions in real arithmetic, but curves introduce interpolation error and function evaluations introduce arithmetic error.

For f with continuous second derivative and |f″|≤M₂ on[a,b], the trapezoidal truncation bound is (b−a)M₂h²/12. One derivation bounds linear interpolation error on each width-h interval by M₂(t−left)(right−t)/2. The interpolation formula follows from repeated Rolle applied to f minus its endpoint line and a multiple of the endpoint-vanishing quadratic. Integrating the bound gives M₂h³/12 per interval, and summing n intervals gives the result. Convex f lies below each chord, so trapezoids overestimate its integral; the sign claim uses convexity, while the magnitude bound uses the derivative control.

Simpson’s composite rule requires positive even n and uses h/3 times endpoint values plus4 times odd interior values and2 times even interior values. Each two-subinterval panel integrates its quadratic interpolant. It is also exact for cubics, since a midpoint-centred cubic’s odd part integrates and samples symmetrically to zero while the even constant/quadratic parts are exact. For continuous fourth derivative bounded by M₄, the standard remainder theorem gives truncation bound (b−a)M₄h⁴/180. The sharp constant is a quadrature remainder result, stated here under its smoothness hypothesis; it is documented in NIST’s quadrature reference.

The fourth-order result cannot be attached to an integrand with an unexamined kink or singularity. Splitting a known kink into separate smooth intervals can restore suitable bounds on each piece; a singular improper endpoint may require a substitution, explicit tail bound or specialised method. Applying an endpoint rule directly to1/sqrt(x) at zero would attempt to evaluate an infinite value, even though its improper integral converges. Mathematical integrability and suitability for a particular sampling formula are different questions.

An error estimate from comparing n and2n values is empirical unless further assumptions link their difference to the true error. A narrow feature can be missed by both grids, just as in Module15. A derivative bound and maximum mesh size can certify unsampled variation, while two matching sums alone cannot. Adaptive methods subdivide where a local estimate indicates difficulty, but they still require an error model and can miss a feature that every inspected sample avoids. Record what information underlies a claimed tolerance.

Arithmetic error adds to truncation. If every sampled value has absolute error at most E, the trapezoidal weights are nonnegative and sum to b−a, so value-error contribution is at most(b−a)E. Simpson weights also sum to b−a for positive even n, giving the same value-error bound. Summation and weight computation can add more error. Refining indefinitely eventually offers diminishing benefit if function evaluation or accumulation dominates, so a convergence table should inspect both trend and plateau.

For a signed integral near zero, relative error can be misleading or undefined because large positive and negative contributions cancel. Use a stated absolute tolerance and, where appropriate, a scale based on ∫|f| or other problem information. Numerical quadrature evaluates an approximation to the mathematical integral; it does not establish an event model, density nonnegativity or differential-equation stability by itself. Those contracts need their own checks.

A first-order ordinary differential equation specifies x′(t)=F(t,x(t)). A solution is a differentiable function on a stated time interval satisfying the equation there; an initial condition x(t₀)=x₀ selects a candidate trajectory. An ODE is not simply a recurrence or a list of samples. Existence and uniqueness depend on conditions on F. This lesson solves selected examples explicitly rather than claiming every right-hand side and initial value has a unique global solution.

For x′=−ax with constant a, direct differentiation verifies x(t)=C exp(−at). The initial condition at zero gives C=x₀. Uniqueness follows without dividing by x: the product exp(at)x(t) has derivative exp(at)[ax+x′]=0, so it is constant on the interval. This includes the equilibrium solution x≡0. For a>0, magnitudes decay and a nonnegative initial state stays nonnegative; a=0 preserves the state, while a<0 causes growth for nonzero initial state. The sign assumption on a determines the dynamical interpretation.

Separation is another route when division is legitimate: x′/x=−a integrates to ln|x|=−at+C on intervals where x≠0. But dividing by x excludes the zero solution, so recover it separately from the original equation. For a general separable x′=g(x)r(t), division by g(x) similarly excludes equilibria at roots of g. Check those roots directly and specify branches and domains before inverting an implicit antiderivative. Transformations can remove valid solutions or introduce implicit expressions outside the original interval.

Euler uses the derivative at the start of each step to make a tangent prediction. With sufficiently smooth exact solution, Taylor gives one-step local discrepancy of order h². Repeated propagation over a fixed final interval can produce global error of order h under appropriate Lipschitz and stability assumptions. For this scalar linear example the formula proves the limiting behaviour directly: with T=nh and h=T/n, the endpoint is x₀(1−aT/n)^n→x₀exp(−aT). The convergence statement holds for fixed T as n→∞, and differs from letting k→∞ at a fixed potentially unsuitable h.

For a>0 and h>0, Euler’s homogeneous factor q=1−ah tends to zero over repeated steps exactly when |q|<1, equivalently0<ah<2. At ah=2 it alternates without shrinking for nonzero initial state. Above two it grows in magnitude with alternating sign. At ah=0 it stays constant. The positive continuous decay is therefore numerically unstable for a sufficiently large step, even though the differential equation itself is stable. The stability region belongs to the method and step, not to the exact model alone.

Nonnegativity is a separate requirement. For a nonnegative initial value, all Euler iterates stay nonnegative when0≤ah≤1; for1<ah<2 they alternate signs while their magnitude decays. Thus a simulation can be asymptotically stable yet violate a physical concentration constraint. At ah=1 the first step goes directly to zero, also a poor finite-time approximation unless that error is acceptable. Check stability, admissible states and accuracy independently instead of using one property as a substitute for the others.

A simulation should report the equation, parameters, initial state, time interval and step rule. If different steps imply different final times, compare each sample with the exact trajectory at its own time. To assess refinement, keep the final time fixed and increase step count. Lab3 does both and explicitly separates long-time stability from finite-time endpoint error. A stable factor cannot make those different experimental questions identical.

These scalar dynamics also prepare for optimisation and AI models. A gradient flow is a continuous equation, while a gradient update is a discrete method with a step-size-dependent stability region. The eigenmode analysis in Module19 extends the scalar factor to quadratic objectives; its assumptions will be stated there. Here the exact solution and explicit recurrence already demonstrate why numerical behaviour must be assessed against the underlying mathematical model.

Lab1 · Quadrature rates and analytic bounds

Download lab1_quadrature.py

"""Composite quadrature on smooth functions with analytic integrals."""
import math

def trapezoid(f, a, b, n):
    if not isinstance(n, int) or n <= 0:
        raise ValueError("n must be a positive integer")
    h = (b-a)/n
    return h*(0.5*(f(a)+f(b))+sum(f(a+i*h) for i in range(1, n)))

def simpson(f, a, b, n):
    if not isinstance(n, int) or n <= 0 or n % 2:
        raise ValueError("n must be a positive even integer")
    h = (b-a)/n
    return h/3*(f(a)+f(b)+sum((4 if i % 2 else 2)*f(a+i*h) for i in range(1, n)))

print("x^4 on [0,1]: exact integral 0.2")
print("n      trapezoid error      Simpson error      trap bound h^2      Simpson bound")
previous = None
for n in [10, 20, 40, 80]:
    h = 1/n
    trap_error = abs(trapezoid(lambda x: x**4, 0, 1, n)-0.2)
    simp_error = abs(simpson(lambda x: x**4, 0, 1, n)-0.2)
    trap_bound, simp_bound = h*h, 24*h**4/180
    print(f"{n:3d} {trap_error:20.10e} {simp_error:18.10e} {trap_bound:18.10e} {simp_bound:18.10e}")
    assert trap_error <= trap_bound+1e-14 and simp_error <= simp_bound+1e-14
    if previous:
        assert trap_error < previous[0] and simp_error < previous[1]
    previous = trap_error, simp_error
print("sin on [0,pi], n=100 Simpson:", simpson(math.sin, 0, math.pi, 100))
assert abs(simpson(math.sin, 0, math.pi, 100)-2) < 2e-8
try:
    simpson(math.sin, 0, math.pi, 9)
except ValueError as error:
    print("Rejected invalid mesh:", error)

x^4 on [0,1]: exact integral 0.2
n      trapezoid error      Simpson error      trap bound h^2      Simpson bound
 10     3.3300000000e-03   1.3333333333e-05   1.0000000000e-02   1.3333333333e-05
 20     8.3312500000e-04   8.3333333334e-07   2.5000000000e-03   8.3333333333e-07
 40     2.0832031250e-04   5.2083333396e-08   6.2500000000e-04   5.2083333333e-08
 80     5.2082519531e-05   3.2552083806e-09   1.5625000000e-04   3.2552083333e-09
sin on [0,pi], n=100 Simpson: 2.0000000108245044
Rejected invalid mesh: n must be a positive even integer

Derive both derivative bounds for x⁴ before running. Compare refinement ratios, note the final rounding allowance, and explain why odd Simpson subinterval counts are rejected. Use the sine example as a second smooth integrand rather than assume the polynomial evidence proves all functions.

Lab2 · Normalisation heights and improper limits

Download lab2_density_and_improper.py

"""Normalisation is an integral condition, distinct from density height."""
import math

def trapezoid(f, a, b, n):
    h = (b-a)/n
    return h*(0.5*(f(a)+f(b))+sum(f(a+i*h) for i in range(1, n)))

for coefficient in [-2.0, 2.0, 3.0]:
    mass = trapezoid(lambda x: coefficient*x, 0, 1, 1000)
    # For this known linear family, c>=0 establishes nonnegativity on [0,1].
    valid = coefficient >= 0 and abs(mass-1) < 1e-12
    print("c / integrated mass / normalised nonnegative candidate:", coefficient, mass, valid)
    assert valid == (coefficient == 2)

mean = trapezoid(lambda x: x*2*x, 0, 1, 1000)
interval_mass = trapezoid(lambda x: 2*x, 0.25, 0.75, 1000)
print("p(1)=2 is a height, not a point probability:", 2.0)
print("Integral x*p(x) / exact weighted mean:", mean, 2/3)
print("Mass on [0.25,0.75] / exact:", interval_mass, 0.5)
assert abs(mean-2/3) < 4e-7 and abs(interval_mass-0.5) < 1e-12

print("epsilon    integral 1/x from epsilon to 1    integral 1/sqrt(x)")
for epsilon in [1e-1, 1e-2, 1e-4, 1e-8]:
    divergent = -math.log(epsilon)
    convergent = 2*(1-math.sqrt(epsilon))
    print(f"{epsilon:.0e} {divergent:32.10f} {convergent:24.10f}")
print("Each cutoff integral is finite; only the limit decides improper convergence.")

c / integrated mass / normalised nonnegative candidate: -2.0 -1.0 False
c / integrated mass / normalised nonnegative candidate: 2.0 1.0 True
c / integrated mass / normalised nonnegative candidate: 3.0 1.5000000000000002 False
p(1)=2 is a height, not a point probability: 2.0
Integral x*p(x) / exact weighted mean: 0.6666669999999999 0.6666666666666666
Mass on [0.25,0.75] / exact: 0.5000000000000001 0.5
epsilon    integral 1/x from epsilon to 1    integral 1/sqrt(x)
1e-01                     2.3025850930             1.3675444680
1e-02                     4.6051701860             1.8000000000
1e-04                     9.2103403720             1.9800000000
1e-08                    18.4206807440             1.9998000000
Each cutoff integral is finite; only the limit decides improper convergence.

Explain the all-domain nonnegativity proof for the linear family, calculate its exact total and interval masses, and compare the weighted integral with2/3. The improper examples use exact antiderivative formulas at finite cutoffs; prove their limiting claims separately.

Lab3 · Decay stability and fixed-time refinement

Download lab3_euler_decay.py

"""A stable continuous decay can have an unstable or negative Euler trace."""
import math

rate = 2.0
print("x'=-2x, x(0)=1; exact x(t)=exp(-2t)")
print("h     Euler factor    x_20          exact at t=20h     behaviour")
for h in [0.1, 0.8, 1.0, 1.1]:
    factor = 1-rate*h
    x = 1.0
    for _ in range(20):
        x += h*(-rate*x)
    behaviour = "decays" if abs(factor) < 1 else "neutral oscillation" if abs(factor) == 1 else "unstable"
    print(f"{h:3.1f} {factor:12.4f} {x:13.7f} {math.exp(-rate*20*h):18.8e} {behaviour}")
    assert abs(x-factor**20) < 1e-10
    if abs(factor) < 1:
        assert abs(x) < 1
    elif abs(factor) > 1:
        assert abs(x) > 1
print("h=0.8 first iterate:", 1-rate*0.8, "(negative despite positive continuous solution)")

print("Fixed final time T=1: n, Euler endpoint, exact endpoint, absolute error")
previous_error = float("inf")
for n in [10, 20, 40, 80]:
    h, x = 1/n, 1.0
    for _ in range(n):
        x += h*(-rate*x)
    exact = math.exp(-rate)
    error = abs(x-exact)
    print(n, f"{x:.10f}", f"{exact:.10f}", f"{error:.10f}")
    assert error < previous_error
    previous_error = error
print("Long-time stability, nonnegativity, and finite-time accuracy are separate checks.")

x'=-2x, x(0)=1; exact x(t)=exp(-2t)
h     Euler factor    x_20          exact at t=20h     behaviour
0.1       0.8000     0.0115292     1.83156389e-02 decays
0.8      -0.6000     0.0000366     1.26641655e-14 decays
1.0      -1.0000     1.0000000     4.24835426e-18 neutral oscillation
1.1      -1.2000    38.3375999     7.78113224e-20 unstable
h=0.8 first iterate: -0.6000000000000001 (negative despite positive continuous solution)
Fixed final time T=1: n, Euler endpoint, exact endpoint, absolute error
10 0.1073741824 0.1353352832 0.0279611008
20 0.1215766546 0.1353352832 0.0137586286
40 0.1285121566 0.1353352832 0.0068231267
80 0.1319378054 0.1353352832 0.0033974778
Long-time stability, nonnegativity, and finite-time accuracy are separate checks.

Classify each factor before running, compare the exact value at each trace’s own final time, and identify the stable sign-violating step. In the second table the final time stays one while the mesh shrinks. Explain what each experiment establishes and what it does not.

Exercises1–12 are required. Optional Exercises13–14 add35 minutes beyond the ten-hour core schedule.

Read integration, improper limits and differential-equation selections in the primary MIT OpenCourseWare calculus course. Consult NIST’s quadrature formulas for the stated trapezoidal and Simpson remainder conventions. The accumulation, substitution, normalisation and decay derivations above are original lesson arguments.

Construct a partition sum, explain why continuity gives integrability, and derive both fundamental-theorem connections. Evaluate a substitution and integration-by-parts example with boundaries. Analyse a convergent and divergent improper limit, normalise a density, and state the quadrature smoothness contracts.

Exit task: x′=−3x,x(0)=2 has exact solution2exp(−3t). Euler is stable for0<h<2/3 and preserves nonnegative states for0<h≤1/3. Compare a stable sign-alternating step with fixed-final-time refinement, and explain why those checks answer different questions.

Ready to move on: you can distinguish continuous mathematical models, accumulated quantities and finite numerical approximations. Multivariable derivatives next extend local sensitivity and chain dependencies to vectors, matrices and automatic differentiation. See the course overview.