Mathematical Model for eliminating pollutes from the great lakes. (I need coding help)

%Matrix that models the flow of pollution through the Great Lakes

a=[(-15/3132), 0, 0, 0, 0; 0, (-38/1322), 0, 0, 0; (15/3132), (38/1322), (-68/782), 0, 0; 0, 0, (68/782),(-85/102), 0; 0, 0, 0, (85/102), (-99/424)]

% The eigan values of this matrix are shown by
eigvalue= eig(a)
   %To find the eigan vectors, set the equation below into MATLAB, where
   %the new matrix "v" shows the eigan vectors(first column is eigenvector to first column eigenvalue)
   [v a]=eig(a);
   v 

     %when m= initial values
     initialval=[3132;1322;782;102;424]
     constants=(v^-1*initialval)% this answer will represent the corresponding constants. All numbers are multiplied by 1000.

So, for a project, we had to develop a mathematical model to drain the pollutants from the great lakes. We put our system of linear equations into a matrix and computed the eigenvalues, eigenvectors, and constants.

Each of the rows represent a flow rate of Pollutants into and out of each lake (positive being into the lake and negative being out of the lake). The rows go down in order of Lake superior, Michigan, Huron, Erie, and Ontario.

My issue is that I have to have a scalar equation that is a function of time but I'm not sure which eigenvalues and eigenvectors correspond to what lakes. Is there a more organized way to do this?

If this isn't enough information, I'd be happy to elaborate. Thank you

Answer 1

[Star Strider]

I don’t know anything about Great Lakes hydrodynamics, and I’m not certain what you are asking or what you want, so this is a guess.

I would approach this as a linear compartmental analysis problem, and use the matrix exponential expm function in a loop.

Example —

a=[(-15/3132), 0, 0, 0, 0; 0, (-38/1322), 0, 0, 0; (15/3132), (38/1322), (-68/782), 0, 0; 0, 0, (68/782),(-85/102), 0; 0, 0, 0, (85/102), (-99/424)];
ic = rand(5, 1);                                    % Initial Concentrations
N = 25;                                             % Iteration Limit
t = linspace(0, 10, N);                             % Time Vector
P = zeros(5, N);                                    % Preallocate

for k1 = 1:N
    P(:,k1) = eye(5) * expm(a*t(k1)) * ic;
end

figure(1)
plot(t, P)
legend('Superior',' Michigan', 'Huron', 'Erie', 'Ontario')
grid

I don’t guarantee that this bears any relation to reality, or to the sort of modeling you need to do.

[Matthew Orsie]

I rewrote the question in hopes that it makes more sense. our resulting equations should look something like (eigenvector)*(constant)e^(eigenvalue*t). My problem is that I don't know how to identify the right values for each lake. I have minimal experience with MATLAB which just kinda spits out numbers and it's up to me to figure out what they are. If there's a way that I can Individually find each eigenvalue and then see it's resulting eigenvector, that would be golden

[Star Strider]

You seem to be describing the Cayley-Hamilton Theorem that is used (among other things) to calculate powers of square matrices, exponentials of square matrices, and related square matrix functions.

My code — using expm — implements the necessary application of the Cayley-Hamilton Theorem and does this for you.

You have a square matrix, so I do not understand the ‘scalar equation’ requirement.

As for the lakes, the rows and columns of your ‘a’ matrix probably sort these for you.

You can implement your own version of expm using the Cayley-Hamilton Theorem (any good linear control textbook, such as T-S Chen’s excellent text, goes into this in requisite detail), use expm, or code your ‘a’ matrix as a system of ordinary differential equations and use ode45 to solve it. I’m lazy, so I just use expm for simple linear problems.