I was wondering if somebody could help me about the following code which doesn't have output as a figure. Thanks in advance.
clcclearn = 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)*x3g=@(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*x1for 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-1f1 = 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-1f1 = 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 figureplot(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')
Best Answer