Trying to write an ODE solver using Backward Euler with Newton-Raphson method

Hi, I'm trying to write a function to solve ODEs using the backward euler method, but after the first y value all of the next ones are the same, so I assume something is wrong with the loop where I use NewtonRoot, a root finding function I wrote previously. Do I need to include a separate loop over the Newton-Raphson method? If so, I'm not sure what index I would use.

function [t,y]=BackwardEuler(F,a,b,y0,N,err,imax)
h=(b-a)/N;
y(1)=y0;
t=a:h:b;
syms f(u)
f(u)=u-y0-h*F(u);
df=diff(f,u);
for ii=1:N
  y(ii+1)=NewtonRoot(f,df,y(ii),err,imax);
end
end

Answer 1

[Abraham Boayue]

Alternatively, you can write a function.

function [t,yb,h] = backwardEuler(a,b,N,M,y0,f,fx)
h   =  (b-a)/N;
t     =  zeros(1,N);
yb    =  zeros(1,N);
t(1)  =  a;
yb(1) =  y0;
for j=1:N
x = yb(j);
t(j+1)=t(j)+h;
for i=1:M
x = x-(x-yb(j)-h*f(t(j+1),x))/(1-h*fx(t(j+1),x));
end
yb(j+1)= x;
end
end

f =  @(t,x) -0.800*x.^(3/2)+10.*2000 *(1 - exp(-3*t));
fx = @(t,x) -0.800*1.5*x.^(1/2);
a = 0;
b = 2;
N = 250;
M = 20;
err = 0.001;
y0 = 2000;
[t,yb,h] = backwardEuler(a,b,N,M,y0,f,fx);
figure()
plot(t,yb,'linewidth',1.5,'color','b')
grid;
a = title('Backward Euler Method');
set(a,'fontsize',14);
a = ylabel('n');
set(a,'Fontsize',14);
a = xlabel('t(s)');
set(a,'Fontsize',14);
axis([0 0.5 0 2000]),

[Abigail Grein]

Do you think there's any way I can use the function I already wrote for Newton's method inside the backward Euler function?

Answer 2

[Abraham Boayue]

I would recommend you writing another function that does the differentiation separately. It would even be best if the function to be differentiated can be done by hand. Can you post the NewtonRoot function? It may give a clue of what's wrong.

[Abigail Grein]

if true
function Xs=NewtonRoot(Fun,FunDer,Xest,Err,imax)
  for i=1:imax
  Xi=Xest-Fun(Xest)/FunDer(Xest); %Xi is the x intercept of the tangent line at Fun(Xest)
  if abs(Fun(Xi))<Err %tolerance in F
      Xs=Xi;
      break
end

[Abigail Grein]

I did think something might have come from trying to make MATLAB do the differentiation I'll try changing that.

[Abigail Grein]

Hmm differentiating by hand didn't change anything. I think it must be the way I'm using the for loop over the NewtonRoot function.

Answer 3

[Abraham Boayue]

Here are two methods that you can use to code Euler backward formula.

%%Method 1
clear variables
close all
a = 0; b = 0.5; 
h = 0.002;
N = (b - a)/h;
M = 20;
n = zeros(1,N);
t = n;
n(1) = 2000;
t(1) = a;
for i=1:N
t(i + 1) = t(i) + h;
x = n(i);
% Newton's implicit method starts.
for j = 1:M
num = x + 0.800*x.^(3/2) *h - 10.*n (1) * (1 - exp(-3*t(i + 1))) *h - n(i) ;
denom = 1 + 0.800*1.5*x.^(1/2) *h;
xnew = x - num/ denom;
if abs((xnew - x)/x) < 0.0001
    break
else
    x = xnew;
end
end
%Newton's method ends.
n(i + 1) = xnew;
end
figure(1)
plot(t,n,'linewidth',1.5,'color','b')
grid;
a = title('Backward Euler Method');
set(a,'fontsize',14);
a = ylabel('n');
set(a,'Fontsize',14);
a = xlabel('t(s)');
set(a,'Fontsize',14);
axis([0 0.5 0 2000]),

%%Method 2
f  =  @(t,x) -0.800*x.^(3/2)+10.*2000 *(1 - exp(-3*t));
fx = @(t,x) -0.800*3/2*x.^(1/2);
a = 0;
b = 2;
n = 250; 
h = (b-a)/n;
y0 = 2000;
t = zeros(1,n);
yb  =zeros(1,n);
t(1)=0;
yb(1)= y0;
% backward Euler method.
for j=1:n
z = yb(j);
t(j+1)=t(j)+h;
for i = 1:20
z = z-(z-yb(j)-h*f(t(j+1),z))/(1-h*fx(t(j+1),z));
end
yb(j+1)=z;
end
figure(2)
plot(t,yb,'linewidth',1.5,'color','r')
grid;
a = title('Backward Euler Method');
set(a,'fontsize',14);
a = ylabel('n');
set(a,'Fontsize',14);
a = xlabel('t(s)');
set(a,'Fontsize',14);
axis([0 0.5 0 2000]),