how to rid of this warning and get correct solution?

2 Ansichten (letzte 30 Tage)
SAHIL SAHOO
SAHIL SAHOO am 14 Okt. 2022
Beantwortet: Gokul Nath S J am 18 Okt. 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

Akzeptierte Antwort

Gokul Nath S J
Gokul Nath S J am 18 Okt. 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.

Weitere Antworten (0)

Community Treasure Hunt

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

Start Hunting!

Translated by