error in ode45. I have declared the vector to be 3x1 but it only reads it as 2x1.
2 visualizzazioni (ultimi 30 giorni)
Mostra commenti meno recenti
Ayush Ranjan
il 21 Apr 2023
Modificato: Stephen23
il 21 Apr 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 Commenti
Risposta accettata
Dyuman Joshi
il 21 Apr 2023
Modificato: Dyuman Joshi
il 21 Apr 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 Commento
Stephen23
il 21 Apr 2023
Modificato: Stephen23
il 21 Apr 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:
Più risposte (0)
Vedere anche
Categorie
Scopri di più su Ordinary Differential Equations in Help Center e File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!