error in ode45. I have declared the vector to be 3x1 but it only reads it as 2x1.
2 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Ayush Ranjan
el 21 de Abr. de 2023
for kpm=27:1:30
tspan = 0:30;
yo=[0.01 0.05 0.179];
[t,y]=ode45(@newwaykumar2004,tspan,yo);
end
function dydt=newwaykumar2004(t,y)
kp=3;
global kpm;
kmp = 25;
klm = 15;
kl = 1;
theta=1;
w=0.5;
function qz=f(m)
qz=1+tanh((m-theta)/w);
end
dpdt=kp*y(1)*(1-y(1))-kpm*y(1)*y(2);
dmdt=(kmp*y(1)+y(3))*y(2)*(1-y(2))-y(2);
dIdt=klm*f(y(2))-kl*y(3);
dydt=[dpdt;dmdt;dIdt];
end
0 comentarios
Respuesta aceptada
Dyuman Joshi
el 21 de Abr. de 2023
Editada: Dyuman Joshi
el 21 de Abr. de 2023
You get the error because dpdt is empty, and thus dydt is 2x1.
Instead of defining kpm as global, it is better to define it as an input to the ode function. MathWorks recommends to avoid the use of global as well
Note that in your code t and y are being overwritten in each iteration of the loop. I have modified t and y to store the result of each iteration
%Defining constant value variables outside the loop
kpm = 27:1:30;
tspan = 0:30;
yo=[0.01 0.05 0.179];
n = numel(kpm);
nt = numel(tspan);
nyo = numel(yo);
%Pre-allocation
t = zeros(nt,n);
y = zeros(nt,nyo,n);
%Ode Loop
for k = 1:n
[t(:,k),y(:,:,k)]=ode45(@(t,y) newwaykumar2004(t,y,kpm(k)), tspan, yo);
end
%t(:,k) represents a solution corresponding to kpm(k)
t
%y(:,:,k) represents a solution corresponding to kpm(k)
y
function dydt=newwaykumar2004(t,y,kpm)
kp=3;
kmp = 25;
klm = 15;
kl = 1;
theta=1;
w=0.5;
%function handle
f = @(m) 1+tanh((m-theta)/w);
dpdt=kp*y(1)*(1-y(1))-kpm*y(1)*y(2);
dmdt=(kmp*y(1)+y(3))*y(2)*(1-y(2))-y(2);
dIdt=klm*f(y(2))-kl*y(3);
dydt=[dpdt;dmdt;dIdt];
end
1 comentario
Stephen23
el 21 de Abr. de 2023
Editada: Stephen23
el 21 de Abr. de 2023
Ayush Ranjan: you should follow Dyuman Joshi's excellent advice to make KPM an input argument (and avoid GLOBAL variables). Note that the ODE45 documentation includes an example showing how to pass extra parameters to the ODE function:
which links in turn to this page:
Más respuestas (0)
Ver también
Categorías
Más información sobre Ordinary Differential Equations en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!