如何使用内置函数 numjac 将 Newton-Raphson 方法应用于 Backward Euler 方法？

Question

我正在编写一个程序，通过 theta 方法求解微分方程组的初值问题。

我的代码如下：

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

根据 https://de.mathworks.com/matlabcentral/newsreader/view_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) 结合使用的逻辑。第一个遵循实数算术序列，而第二个期望整数。
处理这样一个事实，即通常不可能以长度 h 的步长从 a 精确到 b。如果 N = floor((b-a)/h)，那么可能需要最后一步长度 b-a-N*h。或者重新定义h=(b-a)/N.
此时您的代码仅在 S 恰好有 2 个元素且 y 是标量时才有效。

如何使用内置函数 numjac 将 Newton-Raphson 方法应用于 Backward Euler 方法？

How to apply Newton-Raphson method to Backward Euler method using built-in function numjac?

matlab

newtons-method

numerical-methods