problem with an M file, taking a Matrix to a power.

% Filename: Project3PartA
% M-file to
%      - find the population a time [1,2,3,4,5,10,25,50,75]
T = [ 0.00 0.00 0.42 
      0.60 0.00 0.00
      0.00 0.75 0.95 ];
x0 = [ 42 ; 0 ; 95 ]; %Initial Population
x1 = (T^1)*x0
x2 = (T^2)*x0
x3 = (T^3)*x0
x4 = (T^4)*x0
x5 = (T^5)*x0
x10 = (T^10)*x0
x20 = (T^20)*x0
x25 = (T^25)*x0
x50 = (T^50)*x0
x75 = (T^75)*x0

here is my M file. I'm using a Leslie matrix (T) to find a population values at the various times [1,2,3,4,5,10....] but for some reason the file only works for T^(1-20). T^25 and up don't provide the right values. I get

x25 =
   1.0e+03 *
    0.3966
    0.2154
    1.0433

x50 =
   1.0e+04 *
    0.4795
    0.2604
    1.2614

x75 =  
   1.0e+05 *

      0.5797
      0.3148
      1.5250

any idea why its giving the wrong values at 25 and above? or how to fix it?

[Image Analyst]

I fixed your code formatting for you this time, but for next time, read this.

Answer 1

[Image Analyst]

I'm not sure of your formula, but are you sure you want to do a matrix exponeniation instead of an element-by-element exponentiation? Maybe try using .^ instead of ^.

x1 = (T.^1)*x0
x2 = (T.^2)*x0
x3 = (T.^3)*x0
x4 = (T.^4)*x0
x5 = (T.^5)*x0
x10 = (T.^10)*x0
x20 = (T.^20)*x0
x25 = (T.^25)*x0
x50 = (T.^50)*x0
x75 = (T.^75)*x0

[Al Thr]

im sure the formula is right but changing it to .^ doesnt change anything. for values below 25 i get my answer, but for some reason at 25 and up it stops outputting what it should.

[Image Analyst]

That is incorrect. It CERTAINLY does change things. Just look:

T = [ 0.00 0.00 0.42
  0.60 0.00 0.00
  0.00 0.75 0.95 ];
T2m = T * T
T2ebye = T .* T

T2m =
      0                     0.315                     0.399
      0                         0                     0.252
     0.45                    0.7125                    0.9025
T2ebye =
     0                         0                    0.1764
     0.36                      0                     0
     0                    0.5625                    0.9025

Those two matrices are certainly NOT the same. So I'm not sure what you did but I don't think you tried what I said.

[Al Thr]

% Filename: Project3PartA % M-file to

% - find the population a time [1,2,3,4,5,10,25,50,75]

T = [ 0.00 0.00 0.42 0.60 0.00 0.00 0.00 0.75 0.95 ];

x0 = [ 42 ; 0 ; 95 ]; %Initial Population

x1 = (T.^1)*x0

x2 = (T.^2)*x0

x3 = (T.^3)*x0

x4 = (T.^4)*x0

x5 = (T.^5)*x0

x10 = (T.^10)*x0

x20 = (T.^20)*x0

x25 = (T.^25)*x0

x50 = (T.^50)*x0

x75 = (T.^75)*x0

this change, using .^ instead of ^ give the exact same results as my first set of code. and both versions give the same incorrect answers at values above 20

[Image Analyst]

It could be that when you raise a number less than 1 to the 20th power, it becomes so small that it's essentially 0.

[Al Thr]

hm thanks for trying, im just going to hand in the values as presented by Matlab. hopefully they'll accept them.

Answer 2

[Star Strider]

I’m not familiar with your model, but my guess is that it began as a matrix differential equation (here I = eye(size(T))):

xdot = T*x, x(0)  —> matrix Laplace transform —> (s*I - T)*X = x0 —> X = inv(s*I - T)*x0 —> matrix inverse Laplace transform —> x = e^(T*t)*x0

and the code to calculate it becomes:

T = [ 0.00 0.00 0.42 
      0.60 0.00 0.00
      0.00 0.75 0.95 ];
x0 = [ 42 ; 0 ; 95 ]; %Initial Population
t = [1,2,3,4,5,10,25,50,75];
for k1 = 1:length(t)
    x(:,k1) = expm(T*t(k1))*x0;
end

figure(1)
semilogy(t, x)
grid
xlabel('t')
ylabel('x(t)')

You must use the matrix exponential function expm to calculate e^(T*t).

[Al Thr]

ok i think i realized my problem, my code did work but for some reason matlab is giving high values in some weird form of a number*e + a value. i assume this is some sort of weird scientific notation i dont know how to read. for example what is 1.525e+.5?

[Star Strider]

Actually, it is 1.525e+05. This is exponential notation, and means 1.525*10^5.