IMPORTANT:
This code is courtesy of William Yeoh. Again, there are 3 m-files here separated by
dividers. Save the first one as
homework1.m, the second one as phi_a.m and the third
one as phi_b.m.
The first one will call on the rest.
%
William W. Yeoh
%
MEAM 535
%
Homework #1
%
function
homework1()
tstop = 20;
a =
1;
[A,B]
= meshgrid(-2:.3:2, -2:.3:2); % creating the mesh to
plot the contour
phi_a = 1/4*A.^4 + 2/3*A.^3 - 3/2*A.^2 + 1/2*B.^2; %
equation phi for part (a)
phi_b = a*(1-cos(A)) + 1/2*B.^2; % equation phi for part (b)
[px_a,py_a] = gradient(phi_a,.3,.3); % getting the direction of the
gradient
[px_b,py_b] = gradient(phi_b,.3,.3); % getting the direction of the
gradient
%
plotting part(a) in two separate graphs
figure;contour(A,B,phi_a,20);hold on;quiver(A,B,px_a,py_a); hold off;
figure;contour(A,B,phi_a,20);hold on;
y0 =
[1.9; 0.5]; % 1st initial conditions
[t,y_a1]
= ode45('phi_a',tstop,y0);
plot(y_a1(:,1),y_a1(:,2));
y0 =
[-0.5; 1.5]; % 2nd initial conditions
[t,y_a2]
= ode45('phi_a',tstop,y0);
plot(y_a2(:,1),y_a2(:,2));
y0 =
[-0.5; 0]; % 3rd initial conditions
[t,y_a3]
= ode45('phi_a',tstop,y0);
plot(y_a3(:,1),y_a3(:,2));
hold off;
% plotting
part(b) in two separate graphs
figure;contour(A,B,phi_b,20);hold on;quiver(A,B,px_b,py_b); hold off;
figure;contour(A,B,phi_b,20);hold on;
y0 =
[1.9; 0.5]; % 1st initial conditions
[t,y_a2]
= ode45('phi_b',tstop,y0,[],a);
plot(y_a2(:,1),y_a2(:,2));
y0 =
[1.5; 1.5]; % 2nd initial conditions
[t,y_a2]
= ode45('phi_b',tstop,y0,[],a);
plot(y_a2(:,1),y_a2(:,2));
y0 =
[-1.8; -1]; % 3rd initial conditions
[t,y_a3]
= ode45('phi_b',tstop,y0,[],a);
plot(y_a3(:,1),y_a3(:,2));
hold off;
% Save this one as phi_a.m
%
William W. Yeoh
%
Auxiliary function to homework1.m
% Septembre 16, 2003
function
dydt = phi_a(t,y)
dydt = [-y(1)^3 - 2*y(1)^2 + 3*y(1); -y(2)];
% Save this one as phi_b.m
%
William W. Yeoh
%
Auxiliary function to homework1.m
% Septembre 16, 2003
function
dydt = phi_b(t,y,flag,a)
dydt = [-a*sin(y(1)); -y(2)];