Newton Raphson for Backwards Euler on System of differential equations

I have been working on implementing Backwards Euler for a system of differential equations:

function dydt = odeSystem(t,y)

dydt(1,1) = y(2);

dydt(2,1) = -y(1);

end

I have attempted to implement Backwards euler on this system using Newton-Raphson to converge to the future value of

Below you can find my attempt on the solution.

if Method == 'IEM' | Method == 'All'

for i4 = 1:length(x)-1 % Over Time

% Forward Euler for initial guess

k1 = odeSystem(x(i4),y(:,i4));

y_old = y(:,i4) + k1*h ;

% Newton Raphson to estimate n+1

error = 1;

tolerance = 1e-6;

while error >= tolerance

f = odeSystem(x(i4+1),y_old)

J = [0 -1;

1 0];

y_new = y_old - (inv(J)*f)

error = max(abs(y_new - y_old)) % (local maximum) absolute error

y_old = y_new;

end

y(:,i4+1) = y_new;

% Overwrite with backward euler

k1 = odeSystem(x(i4+1),y(:,i4+1));

y(:,i4+1) = y(:,i4) + k1*h;

end

figure(2),clf(2)

sgtitle('Implicit Methods')

subplot(211)

hold on

plot(x,y(1,:),'r')

plot(x,y(2,:),'b')

axis([0 100 -1 1])

legend('Euler d\theta/dt','Euler d2\theta/dt2')

title(['Implicit Euler with h = ',num2str(h)])

grid on;

hold off

end

The Newton-Raphson method does not seem to coverge at all as it explodes. Is my implementation of this Newton-Raphson method wrong?

Thanks in advance!

Answer 1

[Torsten]

  y0(1,1) = 0;
  y0(2,1) = 1;
  x = 0:0.01:2*pi;
  h = x(2)-x(1);
  y = zeros(2,numel(x));
  y(:,1) = y0;
  for i4 = 1:numel(x)-1
     % Forward Euler for initial guess
     k1 = odeSystem(x(i4),y(:,i4));
     y_guess = y(:,i4) + k1*h ;
     % Newton Raphson to estimate n+1
     error = 1;
     tolerance = 1e-6;
     y_new = y_guess;
     y_old = y(:,i4);
     while error >= tolerance
        f  = y_new - y_old - h*odeSystem(x(i4+1),y_new)
        J = [1 -h;
             h 1];
        y_iter = y_new - J\f;
        error = max(abs(y_iter - y_new))                   % (local maximum) absolute error
        y_new = y_iter;
     end
     y(:,i4+1) = y_new;
  end
  plot(x,y(1,:),x,y(2,:))


function dydt = odeSystem(t,y)
  dydt(1,1) = y(2); 
  dydt(2,1) = -y(1);
end