如何在MATLAB中使用numjac函数？求示例及参数说明

Alright, let's tackle this numjac confusion head-on—since the docs are lacking, I'll walk you through its parameters, real-world usage examples, and how it fits into your trapezoidal rule work for computing partial derivatives with respect to y in MATLAB.

First, let's clarify the typical parameter set (since implementations can vary slightly, but this aligns with the 9+ argument structure you mentioned):
The standard call looks like:

J = numjac(fun, y0, t, args, h, ylow, yhigh, reltol, abstol, ...)

Here's what each argument actually does:

• fun: The handle to your ODE function. It needs to follow the format dy = fun(t, y, args)—takes time t, state vector y, and any extra parameters args, outputs the state derivative dy.
• y0: The specific state point where you want to compute the Jacobian (your current/initial y for the trapezoidal step).
• t: The time corresponding to y0 (critical for time-dependent systems, even if you're focusing on y partial derivatives right now).
• args: Any extra parameters to pass to your ODE function (e.g., coefficients, constants). If you don't have any, pass [].
• h: The step size for numerical differentiation. A safe default is sqrt(eps) (MATLAB's machine epsilon square root), but you can tweak it for precision/speed.
• ylow/yhigh: Bounds for your state variables (if you have constraints like non-negative concentrations). Use -inf and inf vectors if there are no bounds.
• reltol/abstol: Relative and absolute error tolerances to control the numerical precision of the Jacobian. Common values are 1e-6 (relative) and 1e-8 (absolute).

Let's say your ODE is dy/dt = -k*y + t (a simple linear system). First, define your ODE function:

function dy = myODE(t, y, args)
    k = args; % Extract the coefficient from args
    dy = -k*y + t;
end

Now, compute the Jacobian (which is the partial derivative of dy/dt with respect to y) at y0=1, t=0, with k=2:

% Set up your parameters
y0 = 1;
t = 0;
args = 2;
h = sqrt(eps);
ylow = -inf;
yhigh = inf;
reltol = 1e-6;
abstol = 1e-8;

% Call numjac
J = numjac(@myODE, y0, t, args, h, ylow, yhigh, reltol, abstol);

% Check the result: The theoretical partial derivative dy/dy is -2, so J should be close to this
disp('Jacobian (dy/dy):');
disp(J);

If you're working with a 2D system (like a harmonic oscillator: dy1/dt = y2, dy2/dt = -y1), here's how it works:

function dy = my2DODE(t, y, ~)
    dy = [y(2); -y(1)]; % No extra args, so we ignore the third input
end

% Define state and time
y0 = [1; 0]; % Initial state: y1=1, y2=0
t = 0;
args = []; % No extra parameters
h = sqrt(eps);
ylow = [-inf; -inf];
yhigh = [inf; inf];
reltol = 1e-6;
abstol = 1e-8;

% Compute Jacobian
J = numjac(@my2DODE, y0, t, args, h, ylow, yhigh, reltol, abstol);

% Theoretical Jacobian is [[0, 1], [-1, 0]]—your numerical result should match this closely
disp('2D System Jacobian:');
disp(J);

For your trapezoidal rule implementation, the Jacobian J (dy/dy) is used to build the system you need to solve for the next state. The trapezoidal update step often requires solving:

(I - (h_step/2)*J) * delta_y = h_step * (f(t, y) + f(t+h_step, y+delta_y))/2

Where I is the identity matrix, h_step is your trapezoidal time step, and f is your ODE function. The numjac output gives you the J you need here.

If you're skeptical of numjac's output, you can implement a simple manual numerical Jacobian to cross-check:

function J = manual_jac(fun, y0, t, args, h)
    n = length(y0);
    J = zeros(n);
    dy0 = fun(t, y0, args); % Base derivative
    for i = 1:n
        y_pert = y0;
        y_pert(i) = y_pert(i) + h; % Perturb one state variable
        dy_pert = fun(t, y_pert, args);
        J(:,i) = (dy_pert - dy0)/h; % Finite difference
    end
end

Compare this function's output to numjac—they should be nearly identical (within your tolerance values).

内容的提问来源于stack exchange，提问作者OliviaHarmon