Plotting and Evaluating PROPT Results

After solving a problem with PROPT, you will probably want to evaluate and/or plot the results.

In PROPT, state and control variables (tomStates and tomControls) are functions of a set of coefficients, and of time. When numeric values are given for the coefficients and for time, states and controls get numeric values.

Example

For example, let's create a couple of states:

>> toms t
>> setPhase(tomPhase('p',t,0,2*pi,10))
>> tomStates x y

We can see how x is constructed by converting it into m-code:

>> [code, tempD] = mcode(x)
 
code =
 out      = interp1p(tempD{1},x_p,t);
 
tempD = 
   [11x1 double]

So x depends on time t and on the coefficient vector x_p. When solving for x in a set of ODE:s or DAE:s, we are actually solving for the coefficients x_p.

For example, we can solve a simple set of ODE:s like this:

>> solution = ezsolve(0,{collocate({dot(x)==y,dot(y)==-x}),initial({x==1,y==0})})

solution = 
   x_p: [11x1 double]
   y_p: [11x1 double]

(Note that the solution does not contain values for t. That is because t is not an unknown t.)

Once we have a solution, we can substitute that into any symbolic expression. For example:

>> x_sol = subs(x,solution)
>> mcode(x_sol)

ans =
 out      = interp1p(tempD{2},tempD{1},t);

Now, x_sol is a symbolic expression that only depends on t. If we substitute in a numeric value for t, then we get a numeric result:

>> subs(x_sol,t==0.2)

ans =
   0.9801

It is possible to substitute in a vector of values, one at a time, using atPoints

>> atPoints(1:6,x_sol)

ans =
   0.1153
  -0.8754
  -1.0613
  -0.2714
   0.7680
   1.1013

Some common substitutions have shortcuts, for example initial(x_sol) sets t to the inital point on the active tomPhase, and collocate(x_sol) uses the collocation points.

The command collocate(t) returns a list of the collocation points.

In order to plot a tomState, one can of course use atPoints to create a vector of values to plot, but there is also a short-cut, called ezplot. As soon as the solution has been computed we can do:

>> ezplot([x,y])