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

Question

I am working on a program that solves the initial value problem for a system of differential equations via the theta method.

My code is as follows:

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

However, the part where I want to calculate the new value of y using Newton's Method, is not working properly. I think this is because the derivate, newton's method needs, is not inserted properly.

I want to use numjac to calculate this derivate. What do the inputs for numjac need to be?

score 0 · Accepted Answer · answered May 04 '15 at 22:15

According to https://de.mathworks.com/matlabcentral/newsreader/view_thread/23776, numjac is exclusively intended to be used for ODE functions of the format doty = f(t,y). And n_df should be the derivative of n_f, thus use the derivative of f in the iteration point x. Thus

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);

There are other strange things in your code, if the above does not result in working code then:

Check the logic of using i in a:h:b in conjunction with T(i) and Y(i). The first follows an arithmetic sequence of real numbers while the second expects an integer.
Deal with the fact that it often is not possible to go from a to exactly b in steps of length h. If N = floor((b-a)/h), then one might need a last step of length b-a-N*h. Or redefine h=(b-a)/N.
Your code at this time only works if S has exactly 2 elements and y is scalar.

Thanks for your tips, I will make changes in my code! – Peter Lawrence May 05 '15 at 13:27 — Peter Lawrence, May 05 '15 at 13:27

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

1 Answers1