ex33ch2.m

these codes are for solving OED with matlab

function ex33ch2
% Define the physical constants:
global m1 m2 L1 L2 g
m1 = 937.5;
m2 = 312.5;
L1 = 16;
L2 = 16;
g = 32;
% First solve the linear model:
options = odeset('Mass',@massL,'MassSingular','no','Events',@events);
[tL,yL,teL,yeL,ieL] = ode45(@odesL,[0 2*pi],[-0.5; -1; 1; 2],options);
if ~isempty(ieL)
    fprintf('Linear model:    lower pallet lost at t = %6.4f.\n',...
             teL(end));
end
plot(tL,yL(:,1),'o',tL,-0.5*cos(2*tL) - 0.5*sin(2*tL))
legend('Computed','Analytical')
title('Linear Model: \theta_1(t)')
figure
plot(tL,yL(:,3),'o',tL,cos(2*tL) + sin(2*tL))
legend('Computed','Analytical',2)
title('Linear Model: \theta_2(t)')
% Now solve the nonlinear model:
options = odeset('Mass',@massN,'MassSingular','no','Events',@events);
[tN,yN,teN,yeN,ieN] = ode45(@odesN,[0 2*pi],[-0.5; -1; 1; 2],options);
if ~isempty(ieN)
    fprintf('Nonlinear model: lower pallet lost at t = %6.4f.\n',...
             teN(end));
end
figure
plot(tL,yL(:,1),'b',tN,yN(:,1),'r')
legend('Linear model','Nonlinear model')
title('Computed \theta_1(t)')
figure
plot(tL,yL(:,3),'b',tN,yN(:,3),'r')
legend('Linear model','Nonlinear model',2)
title('Computed \theta_2(t)')
%==================================================================
function dydt = odesL(t,y)
% Linear model:
global m1 m2 L1 L2 g
dydt = zeros(4,1);
dydt(1) = y(2);
dydt(2) = -(m1 + m2)*g*y(1);
dydt(3) = y(4);
dydt(4) = -m2*g*y(3);

function A = massL(t,y)
global m1 m2 L1 L2 g
A = zeros(4);
A(1,1) = 1;
A(2,2) = (m1 + m2)*L1;
A(2,4) = m2*L2;
A(3,3) = 1;
A(4,2) = m2*L1;
A(4,4) = m2*L2;

function dydt = odesN(t,y)
% Nonlinear model:
global m1 m2 L1 L2 g
dydt = zeros(4,1);
dydt(1) = y(2);
dydt(2) = m2*L2*y(4)*y(4)*sin(y(3) - y(1)) - (m1 + m2)*g*sin(y(1));
dydt(3) = y(4);
dydt(4) = -m2*L1*y(2)*y(2)*sin(y(3) - y(1)) - m2*g*sin(y(3));

function A = massN(t,y)
global m1 m2 L1 L2 g
A = zeros(4);
A(1,1) = 1;
A(2,2) = (m1 + m2)*L1;
A(2,4) = m2*L2*cos(y(3) - y(1));
A(3,3) = 1;
A(4,2) = m2*L1*cos(y(3) - y(1));
A(4,4) = m2*L2;

%==================================================================
function [value,isterminal,direction] = events(t,y)
global m1 m2 L1 L2 g
value = y(3) - pi/3;
isterminal = 1;
direction  = 0;