Trying to write an ODE solver using Backward Euler with Newton-Raphson method
29 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Abigail Grein
el 13 de Abr. de 2018
Comentada: Abigail Grein
el 15 de Abr. de 2018
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
0 comentarios
Respuesta aceptada
Abraham Boayue
el 15 de Abr. de 2018
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]),
Más respuestas (2)
Abraham Boayue
el 13 de Abr. de 2018
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.
3 comentarios
Abraham Boayue
el 15 de Abr. de 2018
Editada: Abraham Boayue
el 15 de Abr. de 2018
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]),
0 comentarios
Ver también
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!