2D ODE with constant? how to solve

9 次查看（过去 30 天）
mays rashad
mays rashad 2021-8-4
回答： Sulaymon Eshkabilov 2021-8-4

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

采纳的回答

Sulaymon Eshkabilov
Ran in:
Most parts of your code is ok, but within the loop, you have overlooked sth and thus, you final solutions are not quite accurate. Here is ODE45 simulation which can be compared with your simulation results.
ICs=[0.6;0.6];
a=0.10;
b=10;
t=[0,60];
F = @(t, z)([a-z(1)+z(1).^2*z(2);b-z(1).^2*z(2)]);
OPTs = odeset('reltol', 1e-6, 'abstol', 1e-9);
[time, z]=ode45(F, t, ICs, OPTs);
figure(2)
plot(time,z(:,1),'b',time,z(:,2),'r')
xlabel('time')
ylabel('x(t) y(t)')
legend('x(t)', 'y(t)', 'location', 'best')
title('Schnackenberg eqn simulation'), xlim([0, 5])
figure(1)
plot(z(:,1),z(:,2),'k')
title('Simulation using ODE45'), grid on
xlabel('x(t)')
ylabel('y(t)')

更多回答（1 个）

Sulaymon Eshkabilov
Use odex (ode23, ode45, ode113, etc.) solvers. See this doc how to employ them in your exercise: https://www.mathworks.com/help/matlab/ref/ode45.html?searchHighlight=ode45&s_tid=srchtitle
  1 个评论
mays rashad
mays rashad 2021-8-4
Is this solution correct?
%x'=a-x+x^2y y'=b-x^2y
clear all,close all, clc
x(1)=0.6;
y(1)=0.6;
a=0.10;
b=10;
h=0.02;
t=0:h:60;
for i=1:(length(t)-1)
k1=h*(a-x(i)+y(i)*x(i)^2);
L1=h*(b-y(i)*x(i)^2);
k2=h*(a-(x(i)+k1/2)+(y(i)+L1/2)*(x(i)^2+k1/2));
L2=h*(b-(y(i)+L1/2)*(x(i)^2+k1/2));
k3=h*(a-(x(i)+k2/2)+(y(i)+L2/2)*(x(i)^2+k2/2));
L3=h*(b-(y(i)+L2/2)*(x(i)^2+k2/2));
k4=h*(a-(x(i)+k3)+(y(i)+L3)*(x(i)^2+k3));
L4=h*(b-(y(i)+L3)*(x(i)^2+k3));
x(i+1)=x(i)+(k1+2*k2+2*k3+k4)*(h/6);
y(i+1)=y(i)+(L1+2*L2+2*L3+L4)*(h/6);
end
plot(t,x,'b',t,y,'r')
xlabel('time')
ylabel('x in blue and y in red')
figure
plot(x,y,'g')
title('2D figure(RK4)')
xlabel('X')
ylabel('Y')

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by