Integration, ODEs and Optimization
|
This content was generated with the assistance of AI and should be verified against the official documentation before being relied on in production. The reference material consulted while preparing these pages is listed in the section bibliography. |
Written against MATLAB R2026b. The sessions in the "Octave 11.3.0 output" blocks were run in GNU Octave 11.3.0 with code that Octave also accepts (integral, trapz, ode45, fzero, fminsearch and the others exist there). Octave has no ode object, so the object-style solver and the automatic-differentiation option are labelled "not run here" and were checked line by line against the documentation. Octave’s solvers choose their own step sequence, so the number of steps they take is not MATLAB’s and no step counts are quoted as MATLAB results.
Three families of numerical tasks live in base MATLAB: integrating a function or a table of samples, solving ordinary differential equations (ODEs) and finding roots and minima. The documentation groups them under "Numerical Integration and Differential Equations" and "Optimization". None of the functions on this page needs a toolbox; the last section says where Optimization Toolbox begins (see Toolboxes and Add-On Products). The running example is the heatflow plate (see The heatflow scenario) seen as a single lumped temperature that cools towards the room.
Integrating a function
integral computes a definite integral of a function handle by adaptive quadrature (see Anonymous Functions and Function Handles).
q = integral(fun,xmin,xmax) q = integral(fun,xmin,xmax,Name=Value)
fun must be vectorized: it receives an array of points and returns an array of values of the same size, so the element-wise operators .*, ./ and .^ are needed. A scalar-only function is the commonest cause of a size error here. The name-value arguments are AbsoluteTolerance (default 1e-10) and RelativeTolerance (default 1e-6), ArrayValued for a function that takes a scalar and returns a fixed-size array, and Waypoints, points the integrator should include in its initial mesh (useful for a kink or a singularity inside the interval). The tolerance options were called AbsTol and RelTol before R2026b; the old names still work. The infinite limits -Inf and Inf are accepted.
Double and triple integrals follow the same pattern:
q = integral2(fun,xmin,xmax,ymin,ymax) q = integral3(fun,xmin,xmax,ymin,ymax,zmin,zmax)
fun takes as many arrays as there are variables, all of the same size. In integral2 the limits ymin and ymax may be function handles of x, which gives non-rectangular regions; the Method option takes auto (default), tiled or iterated.
f = @(x) exp(-x.^2);
q = integral(f, 0, 1);
fprintf('integral = %.10f\n', q);
fprintf('erf check = %.10f\n', sqrt(pi)/2*erf(1));
q2 = integral2(@(x,y) x.*y, 0, 1, 0, 2);
fprintf('integral2 = %.10f\n', q2);
q3 = integral3(@(x,y,z) x+y+z, 0, 1, 0, 1, 0, 1);
fprintf('integral3 = %.10f\n', q3);
% area of the triangle 0 <= y <= x <= 1: the upper limit in y is a handle
qd = integral2(@(x,y) ones(size(x)), 0, 1, 0, @(x) x);
fprintf('triangle = %.10f\n', qd);
qa = integral(@(x) [x, x.^2], 0, 1, 'ArrayValued', true)Octave 11.3.0 output (MATLAB’s layout of the last line may differ):
integral = 0.7468241328
erf check = 0.7468241328
integral2 = 1.0000000000
integral3 = 1.5000000000
triangle = 0.5000000000
qa =
0.5000 0.3333Integrating samples: trapz and cumtrapz
When you have measurements instead of a formula, integrate the samples with the trapezoidal rule.
Q = trapz(Y) Q = trapz(X,Y) Q = trapz(___,dim) Q = cumtrapz(X,Y)
trapz returns the integral over the whole vector (or along dim, by default the first dimension whose size is not 1); cumtrapz returns the running integral, with the same size as Y. X may be a vector of sample coordinates, unevenly spaced if need be, or a scalar spacing. In heatflow, trapz(time, power) turns a logged heater power into energy, and cumtrapz gives the energy delivered so far.
t = [0 1 2 3 4 5];
y = [0 1 4 9 16 25]; % y = t.^2, exact integral over [0,5] is 125/3
fprintf('trapz = %g\n', trapz(t, y));
disp(cumtrapz(t, y))
x = linspace(0,1,101);
fprintf('trapz exp(-x^2), 101 points = %.8f\n', trapz(x, exp(-x.^2)));
A = [1 2; 3 4];
disp(trapz(A, 1)), disp(trapz(A, 2))Octave 11.3.0 output:
trapz = 42.5
0 0.5000 3.0000 9.5000 22.0000 42.5000
trapz exp(-x^2), 101 points = 0.74681800
2 3
1.5000
3.5000The trapezoid result 42.5 against the exact 41.67 shows the rule’s error for a curved function on a coarse grid; the second line shows 101 points reproducing the integral value of the previous section to about five digits. Prefer integral when you can evaluate the function where you like, and trapz when the samples are all you have.
Ordinary differential equations: ode45, ode23 and ode15s
The solvers handle initial value problems of the form , where the unknown is a vector y and you give its value at the start time. A higher-order equation is rewritten as a first-order system, with one state per derivative.
[t,y] = ode45(odefun,tspan,y0) [t,y] = ode45(odefun,tspan,y0,options) [t,y,te,ye,ie] = ode45(odefun,tspan,y0,options) sol = ode45(___)
odefun is a function handle @(t,y) returning the derivative as a column vector. tspan of two elements [t0 tf] lets the solver choose the output times; a longer vector requests output at exactly those times. The same call forms exist for ode23 and ode15s, among others. Which to try first:
| Solver | Use |
|---|---|
|
A Runge-Kutta (4,5) pair. The documentation calls it the first solver to try for most nonstiff problems, at medium accuracy. |
|
A lower-order pair, better for crude tolerances or mildly stiff problems. |
|
Nonstiff solvers for strict tolerances, expensive functions or very smooth solutions. |
|
Variable-order solver for stiff problems and differential-algebraic equations. Try it when |
|
Other stiff, moderately stiff and fully implicit solvers. |
A problem is stiff when the solution has components that change on very different time scales; the documentation notes that stiffness has no precise definition, and that a nonstiff solver which crawls forward with tiny steps is the typical symptom.
Options are set with odeset, which returns a structure that is passed as the fourth argument. The commonly used names are RelTol (default 1e-3), AbsTol (default 1e-6), MaxStep, Events, Mass, Jacobian, NonNegative, OutputFcn and Stats. With a single output, sol = ode45(...) returns a structure that deval evaluates anywhere in the interval.
The plate’s cooling curve as an ODE
Treat the plate as one lumped temperature that loses heat to the room at temperature in proportion to the difference (Newton’s law of cooling):
with the exact solution , which gives something to check against. Take , and (all in the heatflow units, degrees and seconds; the values are illustrative).
k = 0.05; Tamb = 20; T0 = 100;
dTdt = @(t, T) -k*(T - Tamb);
exact = @(t) Tamb + (T0 - Tamb)*exp(-k*t);
[t, T] = ode45(dTdt, [0 60], T0);
fprintf('T(end) = %.6f, exact = %.6f\n', T(end), exact(60));
opts = odeset('RelTol', 1e-8, 'AbsTol', 1e-10);
[t2, T2] = ode45(dTdt, [0 60], T0, opts);
fprintf('tight tolerances: error = %.2e\n', abs(T2(end) - exact(60)));
[t3, T3] = ode45(dTdt, 0:10:60, T0); % output at requested times
disp([t3 T3 exact(t3)])
[t5, T5] = ode15s(dTdt, [0 60], T0);
fprintf('ode15s T(end) = %.4f\n', T5(end));Octave 11.3.0 output (the error and the last digits depend on the solver’s step choices, so MATLAB’s values will differ in the trailing digits):
T(end) = 23.982976, exact = 23.982965
tight tolerances: error = 3.83e-08
0 100.0000 100.0000
10.0000 68.5224 68.5225
20.0000 49.4303 49.4304
30.0000 37.8504 37.8504
40.0000 30.8268 30.8268
50.0000 26.5668 26.5668
60.0000 23.9830 23.9830
ode15s T(end) = 23.9670The default tolerances give about five correct digits, tightening them costs more steps, and ode15s, a stiff solver with a lower accuracy at default tolerances, lands near but not on the answer: the choice of solver is a trade between robustness and accuracy, not a question with one right answer.
Events: when does the plate reach 40 degrees? An event function returns three outputs: the value whose zero crossing is sought, whether to stop (isterminal) and the crossing direction (0 for any). The extra outputs te, ye and ie return the event time, state and index. A root finder gives the same answer from the exact curve, and a second-order system shows how to rewrite it as a first-order one.
k = 0.05; Tamb = 20; T0 = 100;
dTdt = @(t, T) -k*(T - Tamb);
ev = @(t, T) deal(T - 40, true, 0); % stop when T crosses 40
opts = odeset('Events', ev);
[t, T, te, Te, ie] = ode45(dTdt, [0 200], T0, opts);
fprintf('event at t = %.4f, T = %.4f; exact t = %.4f\n', te, Te, log(80/20)/k);
tc = fzero(@(t) Tamb + (T0 - Tamb)*exp(-k*t) - 40, [0 200]);
fprintf('fzero t = %.8f\n', tc);
% y'' = -y as the first-order system [y; y']' = [y'; -y]
f = @(t, y) [y(2); -y(1)];
[t, y] = ode45(f, [0 2*pi], [1; 0]);
fprintf('y(end) = [%.4f %.4f]\n', y(end,1), y(end,2));Octave 11.3.0 output:
event at t = 27.7257, T = 40.0000; exact t = 27.7259
fzero t = 27.72588722
y(end) = [0.9999 -0.0001]The harmonic oscillator should return to its start [1 0] after one period; the small error is the solver’s tolerance, not a bug. The event time is accurate to the solver’s interpolation, which is why a root finder on an exact expression is the more precise check.
The ode object and solve (R2023b)
R2023b introduced an object-based interface. You create an ode object, set properties for the problem, and call solve, which returns an ODEResults object. The documentation example solves a single equation like this (not run here, copied in structure from the reference page):
F = ode;
F.ODEFcn = @(t,y) 2*t;
F.InitialValue = 0;
F.SelectedSolver
sol = solve(F,0,10)
plot(sol.Time,sol.Solution,"-o")The properties used on this page are ODEFcn, InitialTime, InitialValue, Jacobian, AbsoluteTolerance, RelativeTolerance, EventDefinition and Solver. Solver is auto by default, which picks a solver from the problem and the tolerances, and the read-only SelectedSolver shows the choice (the value may change when you change another property). It also accepts stiff, nonstiff or a specific name such as ode45, ode15s, ode113, ode78 and ode89; it covers the SUNDIALS solvers and the delay-equation solvers too. EventDefinition takes an odeEvent object instead of a [value,isterminal,direction] function, and related objects describe a Jacobian (odeJacobian), a mass matrix (odeMassMatrix) and parameter sensitivity (odeSensitivity). solve(F,t0,tf) integrates over an interval, solve(F,t) returns values at the given times, and Refine=N refines the output. The results object has the properties Time, Solution and StopReason.
The functional interface (ode45 and friends with odeset) is still documented and is what most existing code and both books use; nothing on its reference pages marks it as removed. Prefer the object when you want the automatic solver choice, named properties instead of a structure of options, or the extra objects for events, sensitivities and Jacobians. Every example below that must run in Octave uses the functional form.
Automatic-differentiation Jacobians (R2026a). A stiff solver needs the Jacobian . The JacobianMethod property of the ode object, introduced in R2026a, sets how the solver obtains it when you do not supply one: numerical (the default) uses finite differences and autodiff uses automatic differentiation, which the reference page says might be faster for large stiff systems and more accurate for sensitivity analyses. If you set both Jacobian and JacobianMethod, the solver does not compute a Jacobian for the object. Automatic differentiation evaluates derivatives numerically at specific values; it does not build a symbolic expression. Since R2026b it supports complex-valued ODEFcn functions. Not run here:
F = ode;
F.ODEFcn = @(t,T) -0.05*(T - 20);
F.InitialValue = 100;
F.Solver = "ode15s";
F.JacobianMethod = "autodiff";
sol = solve(F,0,60);Finding roots and minima
Base MATLAB has a small, dependable set of scalar solvers.
| Function | What it does |
|---|---|
|
A root of a scalar function of one variable. A scalar |
|
A local minimum of a continuous scalar function on a finite interval, by golden-section search and parabolic interpolation. It does not evaluate the endpoints; if the minimum is at one, it returns an interior point near it. |
|
A local minimum of a function of several variables with no constraints, by the Nelder-Mead simplex method. It uses no derivatives, so it suits non-smooth functions, but it is slow in many variables and is not guaranteed to converge. It works with real variables only. |
|
The vector |
Each takes an options structure from optimset as its last input (for example optimset(TolX=1e-10), with names Display, TolX, TolFun, MaxIter, MaxFunEvals, FunValCheck, OutputFcn, PlotFcns); optimset covers exactly these four solvers. They also return fval, an exitflag (1 means converged for fzero) and an output structure; always check exitflag before trusting x. Every result is a local answer: for a function with several minima the one found depends on the starting point.
f = @(x) (x-2).^2 + 1;
[x, fv, flag] = fminbnd(f, 0, 5);
fprintf('fminbnd x = %.5f f = %.5f flag=%d\n', x, fv, flag);
g = @(p) (1-p(1))^2 + 100*(p(2)-p(1)^2)^2; % Rosenbrock valley
[p, fv] = fminsearch(g, [-1.2 1]);
fprintf('fminsearch p = [%.3f %.3f] f = %.2e\n', p(1), p(2), fv);
o = optimset('TolX', 1e-10, 'TolFun', 1e-10);
p = fminsearch(g, [-1.2 1], o);
fprintf('tight p = [%.6f %.6f]\n', p(1), p(2));
[x, fv, flag] = fzero(@cos, [1 2]);
fprintf('fzero x = %.10f (pi/2 = %.10f) flag=%d\n', x, pi/2, flag);
C = [1 0; 0 1; 1 1]; d = [1; -1; 0.5];
[x, resnorm] = lsqnonneg(C, d);
disp(x'), fprintf('resnorm = %.4f\n', resnorm);
% fit a decay constant to noise-free samples with fminsearch
t = (0:5)'; y = 3*exp(-0.7*t);
res = @(b) sum((b(1)*exp(-b(2)*t) - y).^2);
b = fminsearch(res, [1 1], optimset('TolX',1e-10,'TolFun',1e-12,'MaxFunEvals',2000));
fprintf('fit b = [%.4f %.4f]\n', b(1), b(2));Octave 11.3.0 output:
fminbnd x = 2.00000 f = 1.00000 flag=1
fminsearch p = [1.000 1.000] f = 1.09e-09
tight p = [1.000000 1.000000]
fzero x = 1.5707963268 (pi/2 = 1.5707963268) flag=1
0.7500 0
resnorm = 1.1250
fit b = [3.0000 0.7000]The Rosenbrock function has its minimum at [1 1]; tightening TolX and TolFun is what reaches it to six digits. In lsqnonneg the second component of the unconstrained solution would be negative, so the constraint pins it to zero. The last example fits the amplitude and decay rate of a heatflow-style exponential; for a model that is linear in its parameters, the backslash operator is a better fit method (see Linear Algebra).
The boundary with Optimization Toolbox
The Optimization topic of base MATLAB lists only fminbnd, fminsearch, lsqnonneg, fzero, the option functions optimset and optimget, and the Optimize Live Editor task. Anything with constraints, gradients or large scale is a different product. The fminsearch page itself points to fminunc, the gradient-based unconstrained solver in Optimization Toolbox; constrained problems (fmincon), linear and quadratic programming, nonlinear least squares and integer programming are also there, and global search belongs to Global Optimization Toolbox. Those solvers take options from optimoptions rather than optimset; the optimset page says it does not set the options of every Optimization Toolbox solver. Their functions are not used here (see Toolboxes and Add-On Products). The base-MATLAB way to stay inside the product is fminsearch with a penalty or a transformation for simple bounds, or lsqnonneg for nonnegative linear fits.
Pitfalls
-
Non-vectorized integrands:
integral(@(x) x^2, 0, 1)fails becausefunreceives an array andx^2is a matrix power; writex.^2. -
fzeroneeds a sign change: a double root such as that ofx.^2is invisible to it, and an interval whose ends have the same sign is an error. -
Local answers:
fminsearchandfminbndreturn a local minimum; try several start points. -
Stiff problems: if
ode45takes an enormous number of tiny steps, switch toode15sbefore tightening anything else. -
State vector orientation:
odefunmust return a column vector; a row result is a size error. -
Tolerances are not accuracies:
RelTolandAbsTolbound the local step error, not the final error. -
Events: when an event is terminal the integration stops there, and
t(end)is the event time. -
Data versus function:
trapzon coarse samples gives a coarse answer; resample or useintegralon an interpolant when accuracy matters.
Where older books differ
The Primer by Davis and the Pérez López text were written when quad, quadl and dblquad were the integration functions; the current documentation recommends integral, integral2 and integral3, which have an ArrayValued option and a unified tolerance interface, and quad, quadl and dblquad are no longer the recommended route (the R2026b reference pages for the integral functions do not mention them). The tolerance names AbsTol and RelTol of integral became AbsoluteTolerance and RelativeTolerance in R2026b, and the old names still work. The books also predate the ode object and solve (R2023b) and automatic-differentiation Jacobians (R2026a), so they show only the ode45 and odeset forms, which remain valid. Where they show optimset and fminsearch, the usage is unchanged.