How to evaluate a symbolic Jacobian?

Dear Readers,

How do I evaluate this Jacobian for 50 different combinations of (p,c,R), calculate its eigenvalues and then plot the results against p? I've been stuck with this problem for almost three days.

syms c p R c2 p2

outp = (c2*p2)/(((99*R)/100 - 99/100)*(p2/p)^((2*R)/(R - 1))*(c2/c)^(R/(2*(R - 1))) + 99/100);

infl = (((2*c2*p2)/(35*(((99*R)/100 - 99/100)*(p2/p)^((2*R)/(R - 1))*(c2/c)^(R/(2*(R - 1))) + 99/100)) + 2/175)/(((99*R)/100 - 99/100)*(p2/p)^((2*R)/(R - 1))*(c2/c)^(R/(2*(R - 1))) + 99/100) - (99*p2)/100 + (99*p2^2)/100 + (3*((c2*p2)/(((99*R)/100 - 99/100)*(p2/p)^((2*R)/(R - 1))*(c2/c)^(R/(2*(R - 1))) + 99/100) + 1/5)^(20/7))/35 + 1/4)^(1/2) + 1/2;

J = [jacobian(outp, c2), jacobian(outp, p2) ; jacobian(infl, c2), jacobian(infl, p2)];

%%%%%%%%

beta = 0.99 ; v = 21 ; gamma = 350 ; eps = 1 ; g = 0.2 ; sigma = 1; alpha = 0.7;

c_sol = zeros(1,50);
p = linspace(0.8,1.2,50);
for   i = 1 : numel(p)
pi_actual = p(i);
fun = @(c)  (1 - beta)*pi_actual*(pi_actual-1) - (v/(alpha*gamma))*(c+g).^((1+eps)/alpha) - ((1-v)/gamma)*(c+g)*c.^(-sigma);
c_sol(i) = fzero(fun,0.5);
end

c = c_sol; c2 = c_sol; p2 = p;
R = p/beta; % interest 

%%%%%%%

eigVals = zeros(2,numel(p));

for k = 1:numel(p)

    %part where I miss code

      eigVals(:,k)

      figure;
      hold on;
      plot(p, eigVals(:,k)', '*');
  end

[Joshua]

What is 'B' in the following line?

eigVals(:,k) = eig(B);

It is not defined in the code and Matlab gives an error.

[Sean]

B should be my evaluated jacobian matrix (J), however thats the part where I miss code.

Thanks

[Joshua]

Oh ok. Your code is a bit complex for me to understand what exactly is happening. From what I understand, you want to define a function 50 different times using different values of c, p and R. Each time you want to take the Jacobian of the function, extract the eigenvalues, and then make 50 different plot the eigenvalues vs p. If so, your code looks mostly correct, but needs put all together. I assume you have tried combining everything. What exact problems/errors are you having?

[Sean]

So the the matrix J is now a 2 x 2 symbolic matrix. The variables c2 and p2 are actually lags of c and p and in the steady state (infinite time horizon) those values are equal. So the first loop gives me the steady state solutions of c, p and R, in total 50 combinations. The next step is to evaluate the symbolic Jacobian at each point and determine the eigenvalues, in total 50 different evaluations. Since the state vector is two dimensional, there will only be two eigenvalues per time. Finally, I want to plot the solutions against p (eigenvalues on y-axis ; p on x-axis). However plotting TWO y-values against ONE x-value might also be a little bit strange, but that is of later concern.

Thanks ! :)

PS. I hope you can help, this is the last step of my thesis and thus of my graduation!

PPS. I'll be back tomorrow, it's late :)

[Joshua]

Ok I think I understand it better now although I will leave the math to you haha. If one graph of all 50 eigenvalue combinations is what you are looking for, you want to change your loop to be something like below.

figure % opens a window once
for k = 1:numel(p)
      J=functionof(p(k),c(k),R(k))
      eigVals(:,k) = eig(J);
      hold on; % tells the figure to wait
      plot(p(k), eigVals(:,k)', '*');
end

Note that I moved figure to before the loop to create just one window. Also,

plot(p(k), eigVals(:,k)', '*');

will plot the two eigen values from iteration k versus the value of p at iteration k. The line J=... is meant to show that you want to take the Jacobian each time through the loop based on the specific p, c, and R you have for iteration k. Then you would simply take the eigenvalue of that Jacobian and plot the eigenvalues.

One last thing. Be careful with defining variables as symbolic variables vs actual decimal numbers. defining p as a symbolic variable and then later as a number could possible mess you up. Good luck!

[Sean]

Thanks!

But the problem arises as what to fill in for 'functionof'. Putting 'J' there as defined as the Jacobian matrix in the previous line does not work, then what alternatives are there left? It is indeed correct I need one graph for 50 different combinations of (p, c, R). Lastly, surely every new matrix will give TWO eigenvalues, how would one plot that against ONE p value?

Kind regards, Sean

[Joshua]

Sean,

Here is how you plot two eigenvalues vs one p value without a loop. You can do it with one plot command as in the comment above, but doing it in two allows you to color them separately and refer to them in a legend if that is relevant.

E=rand(2,10)
p=1:length(E)
plot(p,E(1,:),'b*')
hold on
plot(p,E(2,:),'r*')
hold on
legend('First eigenvalue','Second eigenvalue')

I don't really have a full answer for the 'functionof', but we can try breaking it down.

1. You are running a for loop for 50 iterations.

2. Make changes to your parameters p, c, and R each time. Are you just changing p, or also c and R? Is p going to be a real number with c and R being symbolic? You will have to figure this math out because I don't really understand the problem well enough.

3. Calculate the jacobian each iteration using the new parameters for the iteration. I would recommend calculating it inside the loop each time for simplicity.

4. Use the 'eig' command to calculate the eigenvalues each time.

5. Plot the eigenvalues each time. Multiple ways to do this, but the following code would be the easiest to put within a loop. Unless c and R are also both changing each iteration, in which case a two dimensional plot would not really work.

plot(p(k), eigVals(1,k)', '*b');
hold on
plot(p(k), eigVals(2,k)', '*r');
hold on

Try making sure you are only doing one thing on each line. Do one iteration by hand and write pseudo code for that process line by line. Then, code each step by itself. Then, put it all in a loop and try it for 50 iteration. I think your biggest problem is deciding when to use symbolic variables, and when to use real numbers. It is probably going to take a lot of trial and error, but that learning to code in a nutshell. MATLAB documentation is very helpful. Good luck.