Looping ode45 Increasing Parameter

I need to loop ode45 and after each iteration a parameter (mp) should increase and then store all solutions in y_end based on ye. What seems to be incorrect in this code? It doesn't seem to be increasing mp.
mp = 0.05:0.01:100; % start at 0.05, increment by 0.01 until 100
y_end = zeros(length(mp)); % final values 4 column 9996 rows?
for ii = 1:length(mp);
Opt = odeset('Events', @myEvent);
[t, y,te,ye,ie] = ode45(@(t, y)odefun(t, y, mp), [0 pi], [O P J K], Opt);
y_end(ii) = ye(1); % the iith entry of y_end should be given by the corresponding ye values, where ye stores 4 variables solved each iteration
end
function [value, isterminal, direction] = myEvent(t, y);
value = double((any(y(1,1:end) >= pi))); % stop at theta ~ pi?
isterminal = 1; % stop the process
direction = 0;
end
function dydt = odefun(t, y, mp) % y = [y z y1 z1]
dydt = zeros(4,1); % A vector of zeros
dydt(1) = y(3); % y_dot = y1
dydt(2) = y(4); % z_dot = z1
dydt(3) = ((g*sin(y(1))*( mcw * lsa - m * lcg - mp * lla)+(mp * lla)*(g*sin(y(2))*cos(y(1)-y(2))-ls*dydt(2).^2*sin(y(1)-y(2))-lla*dydt(1).^2*cos(y(1)-y(2))*sin(y(1)-y(2)))))/((mcw * (lsa)^2) + (1/12 * m * (lla + lsa)^2 + m * (lcg)^2) + mp*(lla).^2*(sin(y(1)-y(2)))^2); % second derivative of theta
dydt(4) = ((lla*((dydt(1).^2*sin(y(1)-y(2))-dydt(3)*cos(y(1)-y(2)))-g*sin(y(2)))))/ls; % second derivative of phi
end

2 Kommentare

Jan
Jan am 21 Jun. 2018
It is impolite, if you remove the question after an answer has been given. It reduces the chance that the forum will care about your questions in the future.
Rena Berman
Rena Berman am 21 Jun. 2018
(Answers Dev) Restored edit

Melden Sie sich an, um zu kommentieren.

Antworten (1)

Star Strider
Star Strider am 19 Jun. 2018

1 Stimme

You did not post your constants, so I cannot run your code.
Most likely, you need to subscript ‘mp’ in your ‘odefun’ call in your ode45 call.
Try this:
[t,y,te,ye,ie] = ode45(@(t, y)odefun(t, y, mp(ii)), [0 pi], [O P J K], Opt);
↑ ← SUBSCRIPT

Kategorien

Mehr zu Loops and Conditional Statements finden Sie in Hilfe-Center und File Exchange

Produkte

Version

R2017b

Gefragt:

am 19 Jun. 2018

Kommentiert:

am 21 Jun. 2018

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by