I think the ode45 code that I wrote works, but how do I plot it together with both forward and backward euler methods?

2 次查看(过去 30 天)
% values of a & b
a = 7;
b = 3;
L = 0.01*b;
R = 10*a;
C = (0.0001)/a;
t = [0;((40*L)/R)];
f = 1/t;
V = 0;
omega = 2*pi*f;
%%% First-order equations
%%% dI = J;
%%% d2I =dJ;
%%% dJ =((V1*omega)/L)*cos(omega*t)-(R/L)*J-(I/(L*C));
%%% Note: R,C & L are considered constants
% Boundary Conditions
ua = 0;
ub = 0;
% Solve system with initital conditions through ode45() solver
mu = 1;
[T,Y] = ode45(@(t,v) [v(2);mu*(1-v(1).^2).*v(2)-v(1)],[0 ((40*L)/R)],[ua;ub]);
% h value given
h = 0.001;
N = round(((40*L)/R)/h);
%%% initial conditions for Euler methods are the same
%%% as boundary conditions
%initial conditions for Euler Method
xf(1) = ua;
yf(1) = ub;
tf(1) = 0;
% Forward Euler method
for n = 1:N
tf(n+1)=tf(n)+h;
xf(n+1)=xf(n)+yf(n)*h;
yf(n+1)=yf(n)+(mu*(1-xf(n)^2)*yf(n)-xf(n))*h;
end
% Backward Euler method
xb(1) = ua;
yb(1) = ub;
tb(1) = 0;
% Backward Euler method
for n = 1:N
tb(n+1) = tb(n)+h;
fun = @(u,v)[(u-xb(n))/h - v;(v-yb(n))/h - (mu*(1-u^2)*v-u)];
sol = fsolve(@(x)fun(x(1),x(2)),[xb(n),yb(n)],optimset('Display','none'));
xb(n+1) = sol(1);
yb(n+1) = sol(2);
end
% Plotting
figure
hold on
plot(T,Y(:,1))
plot(tf,xf)
plot(tb,xb)
hold off

回答(1 个)

Torsten
Torsten 2023-4-28
I changed your code as to make it work (see above).

类别

Help CenterFile Exchange 中查找有关 Ordinary Differential Equations 的更多信息

产品


版本

R2023a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by