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

1 次查看（过去 30 天）
Alvin Sim
Alvin Sim 2023-4-28
回答： Torsten 2023-4-28
Ran in:
% 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 Center 和 File 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