Problem #29

The position of a certain spring-mass system satisfies the initial value problem

(a) Plot the solution of this initial value problem.

(b) Plot u versus t and u' versus t on the same axes.

(c) Plot u' versus u in the phase plane (see Problem 28). Identify several corresponding points on the curves in parts (b) and (c). What is the direction of motion on the phase plane?

Part (a)

The solution to this initial value problem is

.

Here is a plot of this function using MATLAB. Remember to use matrix arithmetic '.*' and './' wherever we want term-by-term multiplication or division.

  >> u=inline('16*exp(-t/8).*sin(sqrt(127)*t/8)/sqrt(127)','t');
  >> figure
  >> hold on
  >> fplot(u,[0,50])
  >> xlabel('t')
  >> ylabel('u(t)')
  >> title('Solution to u'''' + u''/4 + 2u = 0, u(0) = 0, u''(0) = 2')
  >>

Part (b)

The derivative of u(t) is

.

Let's look at a plot of both u(t) and u'(t) on the same axes: 

  >> uprime=inline('2*exp(-t/8).*cos(sqrt(127)*t/8)-2*exp(-t/8).*sin(sqrt(127)*t/8)/sqrt(127)','t');
  >> figure
  >> hold on
  >> fplot(u,[0,50])
  >> fplot(uprime,[0,50],'k')
  >> xlabel('t')
  >> ylabel('Blue is u(t), Black is u''(t)')
  >> title('u(t) and u''(t) versus t')
  >>

In part (c), we are asked to compare points that plot with the one we generate here. Let's put a few markers on the graph for later comparison. We'll use stars on the original function u(t) and plus signs on the derivative u'(t).

  >> plot(0,feval(u,0),'rp')            % red star at (0, u(0))
  >> plot(10,feval(u,10),'mp')          % magenta star at (10, u(10))
  >> plot(20,feval(u,20),'kp')          % black star at (20, u(20))
  >> plot(30,feval(u,30),'gp')          % green star at (30, u(30))
  >> plot(40,feval(u,40),'cp')          % aqua star at (40, u(40))
  >> plot(0,feval(uprime,0),'r+')       % red plus at (0, u'(0))
  >> plot(10,feval(uprime,10),'m+')     % magenta plus at (10, u'(10))
  >> plot(20,feval(uprime,20),'k+')     % black plus at (20, u'(10))
  >> plot(30,feval(uprime,30),'g+')     % green plus at (30, u'(30))
  >> plot(40,feval(uprime,40),'c+')     % cyan plus at (40, u'(40))
  >>

Part (c)

MATLAB plots functions by plotting pairs (x, y), so it is easy to plot parametric functions (x(t), y(t)) to make a phase plot.

  >> figure
  >> hold on
  >> plot(u(0:0.01:50),uprime(0:0.01:50))       % Plot (u(t),u'(t)) for t's between 
  >> xlabel('u(t)')                             %    0 and 50 in increments of 0.01
  >> ylabel('u''(t)')
  >> title('Phase Plot of u'' versus u')
  >>

To compare this graph to the one in part (b) and determine the direction of motion on the phase portrait, we need to mark a few points with stars. The colors here correspond to the same colors on the graph in part (b). 

  >> plot(feval(u,0),feval(uprime,0),'rp')       % red star at (u(0), u'(0))
  >> plot(feval(u,10),feval(uprime,10),'mp')     % magenta star at (u(10), u'(10))
  >> plot(feval(u,20),feval(uprime,20),'kp')     % black star at (u(20), u'(20))
  >> plot(feval(u,30),feval(uprime,30),'gp')     % green star at (u(30), u'(30))
  >> plot(feval(u,40),feval(uprime,40),'cp')     % cyan star at (u(40), u'(40))
  >>

We can see by the progression of stars that this spiral is being traversed clockwise, spiraling inward. This can also be deduced from the graph in part (b), since both functions u(t) and u'(t) decay to zero.