I've tried to make a simulation with matlab but it can't work. especially for the function commands that I have described. I'm not good at matlab. can you help me to check the
1 view (last 30 days)
Show older comments
clc;clear all;format long;
[T y]=ode45('ftuber',0 200,[10000 0 0 2000 0 500 140]',10^-7);
[T0 y0]=ode45('fnaik',0 200,[10000 0 0 2000 0 500 140]',10^-7);
[T1 y1]=ode45('fturun',0 200,[10000 0 0 2000 0 500 140]',10^-7);
figure (1);
plot(T,y(:,1),'.-',T0,y0(:,1),'.-',T1,y1(:,1),'.-');
title('Grafik A')
legend('u=50','u=30','u=25')
xlabel('waktu (hari)');
ylabel('MA(t) (sel)');
figure (2);
plot(T,y(:,1),'.-',T0,y0(:,1),'.-',T1,y1(:,1),'.-');
title('Grafik B')
legend('u=50','u=30','u=25')
xlabel('waktu (hari)');
ylabel('MA(t) (sel)');
figure (3);
plot(T,y(:,1),'.-',T0,y0(:,1),'.-',T1,y1(:,1),'.-');
title('Grafik C')
legend('u=50','u=30','u=25')
xlabel('waktu (hari)');
ylabel('MA(t) (sel)');
figure (4);
plot(T,y(:,1),'.-',T0,y0(:,1),'.-',T1,y1(:,1),'.-');
title('Grafik D')
legend('u=50','u=30','u=25')
xlabel('waktu (hari)');
ylabel('MA(t) (sel)');
figure (5);
plot(T,y(:,1),'.-',T0,y0(:,1),'.-',T1,y1(:,1),'.-');
title('Grafik E')
legend('u=50','u=30','u=25')
xlabel('waktu (hari)');
ylabel('MA(t) (sel)');
function fv=fgliom(T,y)
p1=0.0068;
p2=0.012;
p3=0.002;
i1=4.7*10^-8;
i2=4.7*10^-5;
i3=4.7*10^-8;
c1=510;
c2=510;
a1=510;
a2=510;
a3=510;
s1=1.8*10^-2;
s2=1.8*10^-3;
s3=1.8*10^-3;
v=2;
j=50;
b=0.2;
u=0;
fq = sign(y(5))/2 + 0.5;
gdot=p1*y(1)*(1-(y(1)/c1))-s1*y(1)*(y(2)+y(3))-(i1*y(1)*y(5))/a1+y(1);
fgdot = sign(-gdot)/2 + 0.5;
fv=zeros(5,1);
fv(1)= p1*y(1)*(1-(y(1)/c1))-s1*y(1)*(y(2)+y(3))-(i1*y(1)*y(5))/a1+y(1);
fv(2)= p2*y(2)*(1-(y(2)+y(3))/c2)-s2*y(1)*y(2)-u*fq*y(1)-(i2*y(1)*y(5))/a2+y(2);
fv(3)= p3*y(3)*(1-(y(2)+y(3))/c2)-s3*y(1)*y(3)+u*fq*y(1);
fv(4)= v*gdot*-(fgdot/c1)*y(4)-(i3*y(4)*y(5))/a3+y(4);
fv(5)= j-b*y(5);
end
function fv=fnaik(T0,y0)
p1=0.0068;
p2=0.012;
p3=0.002;
i1=4.7*10^-8;
i2=4.7*10^-5;
i3=4.7*10^-8;
c1=510;
c2=510;
a1=510;
a2=510;
a3=510;
s1=1.8*10^-2;
s2=1.8*10^-3;
s3=1.8*10^-3;
v=2;
j=50;
b=0.2;
u=10^-3;
f(y5)= sign(y(5))/2 + 0.5;
gdot=p1*y(1)*(1-(y(1)/c1))-s1*y(1)*(y(2)+y(3))-(i1*y(1)*y(5))/a1+y(1);
fgdot = sign(-gdot)/2 + 0.5;
fv=zeros(5,1);
fv(1)= p1*y0(1)*(1-(y0(1)/c1))-s1*y0(1)*(y0(2)+y0(3))-(i1*y0(1)*y0(5))/a1+y0(1);
fv(2)= p2*y0(2)*(1-(y0(2)+y0(3))/c2)-s2*y0(1)*y0(2)-u*F*[y0(5)*y0(1)-(i2*y0(1)*y0(5))/a2+y0(2);
fv(3)= p3*y0(3)*(1-(y0(2)+y0(3))/c2)-s3*y0(1)*y0(3)+u*F*(y0(5))*y0(1);
fv(4)= v*gdot*F*-(fgdot/c1)*y0(4)-(i3*y0(4)*y0(5))/a3+y0(4);
fv(5)= j-b*y0(5);
end
function fv=fturun(T1,y1)
p1=0.0068;
p2=0.012;
p3=0.002;
i1=4.7*10^-8;
i2=4.7*10^-5;
i3=4.7*10^-8;
c1=510;
c2=510;
a1=510;
a2=510;
a3=510;
s1=1.8*10^-2;
s2=1.8*10^-3;
s3=1.8*10^-3;
v=2;
j=50;
b=0.2;
u=10^-2;
fq = sign(y(5))/2 + 0.5;
gdot=p1*y(1)*(1-(y(1)/c1))-s1*y(1)*(y(2)+y(3))-(i1*y(1)*y(5))/a1+y(1);
fgdot = sign(-gdot)/2 + 0.5;
fv=zeros(5,1);
fv(1)= p1*y1(1)*(1-(y1(1)/c1))-s1*y1(1)*(y1(2)+y1(3))-(i1*y1(1)*y1(5))/a1+y1(1);
fv(2)= p2*y1(2)*(1-(y1(2)+y1(3))/c2)-s2*y1(1)*y1(2)-u*fq*y1(1)-(i2*y1(1)*y1(5))/a2+y1(2);
fv(3)= p3*y1(3)*(1-(y1(2)+y1(3))/c2)-s3*y1(1)*y1(3)+u*fq*y1(1);
fv(4)= v*gdot*F*-(fgdot/c1)*y1(4)-(i3*y1(4)*y1(5))/a3+y1(4);
fv(5)= j-b*y1(5);
end
2 Comments
Torsten
on 28 Jul 2022
Edited: Torsten
on 28 Jul 2022
The call to the ODE solvers is
[T,Y] = ode45(@f,tspan,y0)
You will have to change your calls accordingly, e.g.
[T y]=ode45(@fgliom,[0 200],vector with five elements as initial conditions for the solution variables);
And the line
f(y5)= sign(y(5))/2 + 0.5;
will throw an error - maybe you mean
f(5)= sign(y(5))/2 + 0.5;
Answers (1)
VBBV
on 28 Jul 2022
Edited: VBBV
on 28 Jul 2022
clc;clear all;format long;
y = [1000 0 0 2000 0]';
[T y]=ode45(@fgliom,[0 200],[10000 0 0 2000 0]',10^-7);
[T0 y0]=ode45(@fnaik,[0 200],[10000 0 0 2000 0]',10^-7);
[T1 y1]=ode45(@fturun,[0 200],[10000 0 0 2000 0]',10^-7);
figure (1);
plot(T,y(:,1),'.-',T0,y0(:,1),'.-',T1,y1(:,1),'.-');
title('Grafik A')
legend('u=50','u=30','u=25')
xlabel('waktu (hari)');
ylabel('MA(t) (sel)');
figure (2);
plot(T,y(:,2),'.-',T0,y0(:,2),'.-',T1,y1(:,2),'.-');
title('Grafik B')
legend('u=50','u=30','u=25')
xlabel('waktu (hari)');
ylabel('MA(t) (sel)');
figure (3);
plot(T,y(:,3),'.-',T0,y0(:,3),'.-',T1,y1(:,3),'.-');
title('Grafik C')
legend('u=50','u=30','u=25')
xlabel('waktu (hari)');
ylabel('MA(t) (sel)');
figure (4);
plot(T,y(:,4),'.-',T0,y0(:,4),'.-',T1,y1(:,4),'.-');
title('Grafik D')
legend('u=50','u=30','u=25')
xlabel('waktu (hari)');
ylabel('MA(t) (sel)');
figure (5);
plot(T,y(:,5),'.-',T0,y0(:,5),'.-',T1,y1(:,5),'.-');
title('Grafik E')
legend('u=50','u=30','u=25')
xlabel('waktu (hari)');
ylabel('MA(t) (sel)');
function fv=fgliom(T,y)
y = [1 1 1 1 0]';
y5 = [1 0 0 2 0]';
p1=0.0068;
p2=0.012;
p3=0.002;
i1=4.7*10^-8;
i2=4.7*10^-5;
i3=4.7*10^-8;
c1=510;
c2=510;
a1=510;
a2=510;
a3=510;
s1=1.8*10^-2;
s2=1.8*10^-3;
s3=1.8*10^-3;
v=2;
j=50;
b=0.2;
u=0;
fq = sign(y(5))/2 + 0.5;
gdot=p1*y(1)*(1-(y(1)/c1))-s1*y(1)*(y(2)+y(3))-(i1*y(1)*y(5))/a1+y(1);
fgdot = sign(-gdot)/2 + 0.5;
fv=zeros(5,1);
fv(1)= p1*y(1)*(1-(y(1)/c1))-s1*y(1)*(y(2)+y(3))-(i1*y(1)*y(5))/a1+y(1);
fv(2)= p2*y(2)*(1-(y(2)+y(3))/c2)-s2*y(1)*y(2)-u*fq*y(1)-(i2*y(1)*y(5))/a2+y(2);
fv(3)= p3*y(3)*(1-(y(2)+y(3))/c2)-s3*y(1)*y(3)+u*fq*y(1);
fv(4)= v*gdot*-(fgdot/c1)*y(4)-(i3*y(4)*y(5))/a3+y(4);
fv(5)= j-b*y(5);
end
function fv=fnaik(T0,y0)
y = [1 1 0 0 0]';
y0 = [1 0 0 2 0]';
p1=0.0068;
p2=0.012;
p3=0.002;
i1=4.7*10^-8;
i2=4.7*10^-5;
i3=4.7*10^-8;
c1=510;
c2=510;
a1=510;
a2=510;
a3=510;
s1=1.8*10^-2;
s2=1.8*10^-3;
s3=1.8*10^-3;
v=2;
j=50;
b=0.2;
u=10^-3;
fy5= sign(y(5))/2 + 0.5;
gdot=p1*y(1)*(1-(y(1)/c1))-s1*y(1)*(y(2)+y(3))-(i1*y(1)*y(5))/a1+y(1);
fgdot = sign(-gdot)/2 + 0.5;
fv=zeros(5,1);
fv(1)= p1*y0(1)*(1-(y0(1)/c1))-s1*y0(1)*(y0(2)+y0(3))-(i1*y0(1)*y0(5))/a1+y0(1);
fv(2)= p2*y0(2)*(1-(y0(2)+y0(3))/c2)-s2*y0(1)*y0(2)-u*fy5*y0(5)*y0(1)-(i2*y0(1)*y0(5))/a2+y0(2);
fv(3)= p3*y0(3)*(1-(y0(2)+y0(3))/c2)-s3*y0(1)*y0(3)+u*fy5*(y0(5))*y0(1);
fv(4)= v*gdot*fy5-(fgdot/c1)*y0(4)-(i3*y0(4)*y0(5))/a3+y0(4);
fv(5)= j-b*y0(5);
end
function fv=fturun(T1,y1)
y = [0 0 1 1 0]';
y1 = [1 0 0 2 0]';
p1=0.0068;
p2=0.012;
p3=0.002;
i1=4.7*10^-8;
i2=4.7*10^-5;
i3=4.7*10^-8;
c1=510;
c2=510;
a1=510;
a2=510;
a3=510;
s1=1.8*10^-2;
s2=1.8*10^-3;
s3=1.8*10^-3;
v=2;
j=50;
b=0.2;
u=10^-2;
fq = sign(y(5))/2 + 0.5;
gdot=p1*y(1)*(1-(y(1)/c1))-s1*y(1)*(y(2)+y(3))-(i1*y(1)*y(5))/a1+y(1);
fgdot = sign(-gdot)/2 + 0.5;
fv=zeros(5,1);
fv(1)= p1*y1(1)*(1-(y1(1)/c1))-s1*y1(1)*(y1(2)+y1(3))-(i1*y1(1)*y1(5))/a1+y1(1);
fv(2)= p2*y1(2)*(1-(y1(2)+y1(3))/c2)-s2*y1(1)*y1(2)-u*fq*y1(1)-(i2*y1(1)*y1(5))/a2+y1(2);
fv(3)= p3*y1(3)*(1-(y1(2)+y1(3))/c2)-s3*y1(1)*y1(3)+u*fq*y1(1);
fv(4)= v*gdot*fq-(fgdot/c1)*y1(4)-(i3*y1(4)*y1(5))/a3+y1(4);
fv(5)= j-b*y1(5);
end
You can use different product multipliers for each function as given in defintion
6 Comments
Torsten
on 28 Jul 2022
I mean the lines
function fv=fgliom(T,y)
y = [1 1 1 1 0]';
y5 = [1 0 0 2 0]';
and
function fv=fnaik(T0,y0)
y = [1 1 0 0 0]';
y0 = [1 0 0 2 0]';
and
function fv=fturun(T1,y1)
y = [0 0 1 1 0]';
y1 = [1 0 0 2 0]';
See Also
Categories
Find more on Introduction to Installation and Licensing 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!