Trying to write an ODE solver using Backward Euler with Newton-Raphson method
54 visualizzazioni (ultimi 30 giorni)
Mostra commenti meno recenti
Abigail Grein
il 13 Apr 2018
Commentato: Abigail Grein
il 15 Apr 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 Commenti
Risposta accettata
Abraham Boayue
il 15 Apr 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]),
Più risposte (2)
Abraham Boayue
il 13 Apr 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 Commenti
Abraham Boayue
il 15 Apr 2018
Modificato: Abraham Boayue
il 15 Apr 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 Commenti
Vedere anche
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!