如何利用内建函数numjac将牛顿-拉夫森方法应用于后向欧拉方法？

Question

我正在做一个程序，通过θ方法来解决微分方程组的初值问题。

我的代码如下：

function [T,Y] = ivpSolver(f, S, y0, theta, h)
%INPUT
%f: a function-handle which describes the right-hand side f(t; u) of
%the differential equation
%S: a vector with requested points in time,
%u0: a column vector with initial conditions %u0
%theta: the parameter θ
%h: the step size h.
%OUTPUT
%T: a column vector of times T
%U: a matrix U with in every row the solution at the time corresponding to
%that row in T.

if length(S) == 2
    a = S(1);
    b = S(2);
end

T = zeros((b-a)/h,1);
Y = zeros((b-a)/h,1);
T(1) = t0;
Y(1) = y0;

tol = eps;
maxiter = 10;

for i = a:h:b
    T(i+1) = T(i) + h;

    %Calculation of new Y using Newton's method
    n_f = @(x) x - Y(i) - h*f(T(i+1),x);
    n_df = @(x)numjac(n_f, 0, Y(i) ,n_f(0, Y(i)), eps);
    ynew = newton(n_f,n_df,Y(i),tol,maxiter);

    Y(i+1) = Y(i) + h * ((1-theta)* f(T(i),Y(i))+ theta* f(T(i+1),ynew));
end

end

然而，我想用牛顿法计算y的新值的部分，是不正确的。我认为这是因为牛顿方法所需要的导数没有被正确地插入。

我想用numjac来计算这个导数。numjac的输入需要什么？

Answer 1

根据thread/23776，numjac专门用于doty = f(t,y)格式的ODE函数。n_df应该是n_f的导数，因此在迭代点x中使用f的导数。因此

n_f = @(x) x - Y(i) - h*f(T(i+1),x);
n_df = @(x) 1 - h*numjac(f, T(i+1), x ,f(T(i+1), x), eps);

在您的代码中还有其他一些奇怪的东西，如果上面的代码不产生工作代码，那么：

• 检查将i in a:h:b与T(i)和Y(i)结合使用的逻辑。前者遵循实数的算术序列，而第二个则期望整数。
• 处理这样一个事实，即通常不可能从a到b的长度为h的步骤。如果是N = floor((b-a)/h)，那么可能需要长度为b-a-N*h的最后一步。或者重新定义h=(b-a)/N。
• 此时，您的代码只有在S只有2个元素且y是标量的情况下才能工作。