Trying to write an ODE solver using Backward Euler with Newton-Raphson method
28 ビュー (過去 30 日間)
古いコメントを表示
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 件のコメント
採用された回答
Abraham Boayue
2018 年 4 月 15 日
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]),
その他の回答 (2 件)
Abraham Boayue
2018 年 4 月 13 日
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 件のコメント
Abraham Boayue
2018 年 4 月 15 日
編集済み: Abraham Boayue
2018 年 4 月 15 日
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 件のコメント
参考
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!