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

Jan
Jan 2018년 6월 21일
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 2018년 6월 21일
(Answers Dev) Restored edit

댓글을 달려면 로그인하십시오.

답변 (1개)

Star Strider
Star Strider 2018년 6월 19일

1 개 추천

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

카테고리

도움말 센터File Exchange에서 Loops and Conditional Statements에 대해 자세히 알아보기

제품

릴리스

R2017b

질문:

2018년 6월 19일

댓글:

2018년 6월 21일

Community Treasure Hunt

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

Start Hunting!

Translated by