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

Question

Abigail Grein 2018-4-13

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/394763-trying-to-write-an-ode-solver-using-backward-euler-with-newton-raphson-method

评论： Abigail Grein 2018-4-15

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 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Abraham Boayue 2018-4-15

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/394763-trying-to-write-an-ode-solver-using-backward-euler-with-newton-raphson-method#answer_315193

在 MATLAB Online 中打开

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]),

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Abigail Grein 2018-4-15

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 2018-4-13

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/394763-trying-to-write-an-ode-solver-using-backward-euler-with-newton-raphson-method#answer_315012

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 个评论
显示 1更早的评论隐藏 1更早的评论

Abigail Grein 2018-4-13

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

Abigail Grein 2018-4-13

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 2018-4-15

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/394763-trying-to-write-an-ode-solver-using-backward-euler-with-newton-raphson-method#answer_315192

编辑：Abraham Boayue 2018-4-15

在 MATLAB Online 中打开

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]),