2D ODE with constant? how to solve

2021 8 4
2 个回答
回答已采纳
更新时间：2021 8 4
2 次查看（30 天）

 采纳的回答

Ran in:

0 个投票

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 个）

0 个投票

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

请先登录，再进行评论。

类别

在 帮助中心 和 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