clear all; close all; clc
An undamped pendulum with no driving force is described by
We reduce this to standard matlab form of a system of first order ODEs by letting and . This leads to:
The phase portrait is a plot of a vector field which qualitatively shows how the solutions to these equations will go from a given starting point. here is our definition of the differential equations:
f = @(t,Y) [Y(2); -sin(Y(1))];
To generate the phase portrait, we need to compute the derivatives and at on a grid over the range of values for and we are interested in. We will plot the derivatives as a vector at each (y1, y2) which will show us the initial direction from each point. We will examine the solutions over the range -2 < y1 < 8, and -2 < y2 < 2 for y2, and create a grid of 20x20 points.
y1 = linspace(-2,8,20); y2 = linspace(-2,2,20); % creates two matrices one for all the x-values on the grid, and one for % all the y-values on the grid. Note that x and y are matrices of the same % size and shape, in this case 20 rows and 20 columns [x,y] = meshgrid(y1,y2); size(x) size(y)
ans = 20 20 ans = 20 20
we need to generate two matrices, u and v. u will contain the value of y1' at each x and y position, while v will contain the value of y2' at each x and y position of our grid.
we preallocate the arrays so they have the right size and shape
u = zeros(size(x)); v = zeros(size(x)); % we can use a single loop over each element to compute the derivatives at % each point (y1, y2) t=0; % we want the derivatives at each point at t=0, i.e. the starting time for i = 1:numel(x) Yprime = f(t,[x(i); y(i)]); u(i) = Yprime(1); v(i) = Yprime(2); end
now we use the quiver command to plot our vector field
quiver(x,y,u,v,'r'); figure(gcf) xlabel('y_1') ylabel('y_2') axis tight equal;
Let's plot a few solutions on the vector field. We will consider the solutions where y1(0)=0, and values of y2(0) = [0 0.5 1 1.5 2 2.5], in otherwords we start the pendulum at an angle of zero, with some angular velocity.
hold on for y20 = [0 0.5 1 1.5 2 2.5] [ts,ys] = ode45(f,[0,50],[0;y20]); plot(ys(:,1),ys(:,2)) plot(ys(1,1),ys(1,2),'bo') % starting point plot(ys(end,1),ys(end,2),'ks') % ending point end hold off
what do these figures mean? For starting points near the origin, and small velocities, the pendulum goes into a stable limit cycle. For others, the trajectory appears to fly off into y1 space. Recall that y1 is an angle that has values from to . The y1 data in this case is not wrapped around to be in this range.
'done' % categories: ODEs % tags: math % post_id = 799; %delete this line to force new post;
ans = done