how to rid of this warning and get correct solution?

2 visualizaciones (últimos 30 días)
SAHIL SAHOO
SAHIL SAHOO el 14 de Oct. de 2022
Respondida: Gokul Nath S J el 18 de Oct. de 2022
ti = 0;
tf = 1E-7;
tspan=[ti tf];
k = 5E-3;
h = 10E-2;
y0= [(h)*rand(2,1); ((-3.14).*rand(1,1) + (3.14).*rand(1,1));
(h)*rand(2,1); ((-3.14).*rand(1,1) + (3.14).*rand(1,1));
(h)*rand(2,1); ((-3.14).*rand(1,1) + (3.14).*rand(1,1));
(h)*rand(2,1); ((-3.14).*rand(1,1) + (3.14).*rand(1,1));
(h)*rand(2,1); ((-3.14).*rand(1,1) + (3.14).*rand(1,1));
((-3.14).*rand(5,1) + (3.14).*rand(5,1))];
yita_mn = [
0 1 0 0 1;
1 0 1 0 0;
0 1 0 1 0;
0 0 1 0 1;
1 0 0 1 0;
]*(k);
N = 5;
tp = 1E-12;
[T,Y]= ode45(@(t,y) rate_eq(t,y,yita_mn,N),tspan./tp,y0);
Warning: Failure at t=4.032134e+01. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (1.136868e-13) at time t.
figure(1)
plot(T,(Y(:,16)),'linewidth',0.8);
hold on
for m = 16:20
plot(T,(Y(:,m)),'linewidth',0.8);
end
hold off
grid on
xlabel("time")
ylabel("phase difference")
set(gca,'fontname','times New Roman','fontsize',18,'linewidth',1.8);
function dy = rate_eq(t,y,yita_mn,N,o)
dy = zeros(4*N,1);
dGdt = zeros(N,1);
dAdt = zeros(N,1);
dOdt = zeros(N,1);
P = 1.25;
a = 5;
T = 2000;
Gt = y(1:3:3*N-2);
At = y(2:3:3*N-1);
Ot = y(3:3:3*N-0);
k = 5E-5;
for i = 1:N
dGdt(i) = (P - Gt(i) - (1 + 2.*Gt(i)).*(At(i))^2)./T ;
dAdt(i) = (Gt(i).*(At(i)));
dOdt(i) = -a.*(Gt(i));
for j = 1:N
dAdt(i) = dAdt(i)+yita_mn(i,j).*(At(j))*sin(Ot(j)-Ot(i));
dOdt(i) = dOdt(i)+yita_mn(i,j).*((At(j)/At(i)))*cos(Ot(j)-Ot(i));
end
end
dy(1:3:3*N-2) = dGdt;
dy(2:3:3*N-1) = dAdt;
dy(3:3:3*N-0) = dOdt;
n1 = (1:5)';
n2 = circshift(n1,-1);
n16 = n1 + 15;
n17 = circshift(n16,-1);
n20 = circshift(n16,1);
j2 = 3*(1:5)-1;
j5 = circshift(j2,-1);
j8 = circshift(j2,-2);
j19 = circshift(j2,1);
dy(n16) = -a.*(Gt(n2)-Gt(n1)) + (k).*(y(j2)./y(j5)).*cos(y(n16)) - (k).*(y( j5)./y(j2)).*cos(y(n16)) + (k).*(y(j8)./y(j5)).*cos(y(n17)) - (k).*(y(j19)./y(j2)).*cos(y(n20));
end

Respuesta aceptada

Gokul Nath S J
Gokul Nath S J el 18 de Oct. de 2022
Dear Sahil Sahoo,
I have tried running your code on my machine. By setting the relative and absolute error values to 1e-2, the warning messages can be avoided.
options = odeset('RelTol',1e-2,'AbsTol',1e-2);
The code should be included before calling ode45 (line 22)
However, at the same time, it is essential to check the validity of the eventual answer.
Also, I have found similar cases that might be helpful. Refer to the link below for the same.

Más respuestas (0)

Categorías

Más información sobre Ordinary Differential Equations en Help Center y File Exchange.

Community Treasure Hunt

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

Start Hunting!

Translated by