how to rid of this warning and get correct solution?
1 view (last 30 days)
Show older comments
SAHIL SAHOO
on 14 Oct 2022
Answered: Gokul Nath S J
on 18 Oct 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);
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
0 Comments
Accepted Answer
Gokul Nath S J
on 18 Oct 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.
0 Comments
More Answers (0)
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!