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 vue (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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 commentaires
Torsten
le 28 Juil 2022
Modifié(e) : Torsten
le 28 Juil 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;
Réponses (1)
VBBV
le 28 Juil 2022
Modifié(e) : VBBV
le 28 Juil 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 commentaires
Torsten
le 28 Juil 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]';
Voir également
Catégories
En savoir plus sur Introduction to Installation and Licensing dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!