No output in my problem

2 vues (au cours des 30 derniers jours)
Mojtaba Mohareri
Mojtaba Mohareri le 18 Mai 2020
I was wondering if somebody could help me about the following code which doesn't have output as a figure. Thanks in advance.
clc
clear
n = 1000;
a = 0;
b = 7;
h = (b-a)/n;
x(1) = 0.009;
y(1) = 0.001;
z(1)=0.99;
t(1)=a;
f=@(t,x1,x2,x3) -0.002*x1+0.018*x2-0.10*x1+(1-0.10)*(0.01+(0.99/(1+exp(-2*t+8)))+1)*((x3*x2)/(x1+x2+x3))+0.06*(1-0.05)*x3
g=@(t,x1,x2,x3) -0.018*x2+0.002*x1-0.10*x2+0.05*0.06*x3+0.10*(0.01+0.99/(1+exp(-2*t+8))+1)*((x3*x2)/(x1+x2+x3))
c=@(t,x1,x2,x3) -(0.01+0.99/(1+exp(-2*t+8))+1)*((x3*x2)/(x1+x2+x3))-0.06*x3+0.10*x2+0.10*x1
for i=1:n-1
t(i+1) = a + i*h;
end
disp('Step t x y z');
disp('----------------------------------------------------');
fprintf('0 %2.2f %2.16f %2.16f %2.16f\n', t(1),x(1),y(1),z(1));
for k=1:n-1
f1 = f(t(k),x(k),y(k),z(k));
g1 = g(t(k),x(k),y(k),z(k));
c1 = c(t(k),x(k),y(k),z(k));
f2 = f(t(k)+(h/2), x(k)+(h/2)*f1 , y(k)+(h/2)*g1 ,z(k)+(h/2)*c1);
g2 = g(t(k)+(h/2), x(k)+(h/2)*f1 , y(k)+(h/2)*g1 ,z(k)+(h/2)*c1);
c2 = c(t(k)+(h/2), x(k)+(h/2)*f1 , y(k)+(h/2)*g1 ,z(k)+(h/2)*c1);
f3 = f(t(k)+(h/2), x(k)+(h/2)*f2 , y(k)+(h/2)*g2,z(k)+(h/2)*c2);
g3 = g(t(k)+(h/2), x(k)+(h/2)*f2 , y(k)+(h/2)*g2,z(k)+(h/2)*c2);
c3 = c(t(k)+(h/2), x(k)+(h/2)*f2 , y(k)+(h/2)*g2,z(k)+(h/2)*c2);
f4 = f(t(k)+h , x(k)+h*f3 , y(k)+h*g3,z(k)+h*c3);
g4 = g(t(k)+h , x(k)+h*f3 , y(k)+h*g3,z(k)+h*c3);
c4 = c(t(k)+h , x(k)+h*f3 , y(k)+h*g3,z(k)+h*c3);
x(k+1) = x(k) + (h/6)*(f1 + 2*f2 + 2*f3 + f4);
y(k+1) = y(k) + (h/6)*(g1+ 2*g2 + 2*g3 + g4);
z(k+1) = z(k) + (h/6)*(c1+ 2*c2 + 2*c3 + c4);
fprintf('%2f %2.2f %2.16f %2.16f %2.16f\n', k, t(k+1),x(k+1),y(k+1), z(k+1));
end
n = 1000;
a = 0;
b = 7;
h = (b-a)/n;
p(1) = 0;
q(1) = 0;
r(1) = 0;
t(1) = b;
f=@(t,x1,x2,x3) 0.018*(x1-x2)+0.10*(x1-x3)+((x3-0.10*x1-(1-0.10)*x2)*(0.01+(0.99)/(1+exp(-2*t+8))+1)*z(k)*(x(k)+z(k)))/((x(k)+y(k)+z(k))^2)
g=@(t,x1,x2,x3) 0.002*(x2-x1)+0.10*(x2-x3)-((x3-0.10*x1-(1-0.10)*x2)*(0.01+(0.99)/(1+exp(-2*t+8))+1)*z(k)*y(k))/((x(k)+y(k)+z(k))^2)
c=@(t,x1,x2,x3) -1+(x3-0.05*x1-(1-0.05)*x2)*0.06+((x3-0.10*x1-(1-0.10)*x2)*(0.01+(0.99)/(1+exp(-2*t+8))+1)*y(k)*(x(k)+y(k)))/((x(k)+y(k)+z(k))^2)
for i=1:n-1
t(i+1) = b - i*h;
end
disp('Step t x y z');
disp('----------------------------------------------------');
fprintf('0 %2.2f %2.16f %2.16f %2.16f\n', t(1),p(1),q(1),r(1));
for k=1:n-1
f1 = f(t(k),p(k),q(k),r(k));
g1 = g(t(k),p(k),q(k),r(k));
c1 = c(t(k),p(k),q(k),r(k));
f2 = f(t(k)-(h/2), p(k)-(h/2)*f1 , q(k)-(h/2)*g1 ,r(k)-(h/2)*c1);
g2 = g(t(k)-(h/2), p(k)-(h/2)*f1 , q(k)-(h/2)*g1 ,r(k)-(h/2)*c1);
c2 = c(t(k)-(h/2), p(k)-(h/2)*f1 , q(k)-(h/2)*g1 ,r(k)-(h/2)*c1);
f3 = f(t(k)-(h/2), p(k)-(h/2)*f2 , q(k)-(h/2)*g2, r(k)-(h/2)*c2);
g3 = g(t(k)-(h/2), p(k)-(h/2)*f2 , q(k)-(h/2)*g2, r(k)-(h/2)*c2);
c3 = c(t(k)-(h/2), p(k)-(h/2)*f2 , q(k)-(h/2)*g2, r(k)-(h/2)*c2);
f4 = f(t(k)-h , p(k)-h*f3 , q(k)-h*g3,r(k)-h*c3);
g4 = g(t(k)-h , p(k)-h*f3 , q(k)-h*g3,r(k)-h*c3);
c4 = c(t(k)-h , p(k)-h*f3 , q(k)-h*g3,r(k)-h*c3);
p(k+1) = p(k) - (h/6)*(f1 + 2*f2 + 2*f3 + f4);
q(k+1) = q(k) - (h/6)*(g1+ 2*g2 + 2*g3 + g4);
r(k+1) = r(k) - (h/6)*(c1+ 2*c2 + 2*c3 + c4);
fprintf('%2f %2.2f %2.16f %2.16f %2.16f\n', k, t(k+1),p(k+1),q(k+1), r(k+1));
end
figure
plot(t,min(max(0,((r-(1/10)*p-(1-(1/10))*q).*z.*y)/(2*(1/100)*(x+y+z))),1))
axis([0 7 0 1]);
xlabel('time')
ylabel('Control U2')

Réponse acceptée

Walter Roberson
Walter Roberson le 18 Mai 2020
plot(t,min(max(0,((r-(1/10)*p-(1-(1/10))*q).*z.*y)./(2*(1/100)*(x+y+z))),1))
That is, your / has to be ./
  1 commentaire
Mojtaba Mohareri
Mojtaba Mohareri le 18 Mai 2020
Thank you very much

Connectez-vous pour commenter.

Plus de réponses (0)

Catégories

En savoir plus sur Graphics Object Programming dans Help Center et File Exchange

Tags

Community Treasure Hunt

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

Start Hunting!

Translated by