error in ode45. I have declared the vector to be 3x1 but it only reads it as 2x1.
2 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Ayush Ranjan
am 21 Apr. 2023
Bearbeitet: Stephen23
am 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 Kommentare
Akzeptierte Antwort
Dyuman Joshi
am 21 Apr. 2023
Bearbeitet: Dyuman Joshi
am 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 Kommentar
Stephen23
am 21 Apr. 2023
Bearbeitet: Stephen23
am 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:
Weitere Antworten (0)
Siehe auch
Kategorien
Mehr zu Ordinary Differential Equations finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!