Homework question, involving MATLAB and least squares problem.

Hello,

I'm new here, I have a class this semester that makes us use a CAS for random problems, even though none of us have ever been taught how to use it. I've been messing around with it for the past couple of weeks and I'm slowly getting the hang of it, but then this week she gives a problem that I'm pretty much lost on what to do with and cannot find help in my notes for the class or the book we use for the class(its a computational physics class).

So I was wondering if any of you could help me figure out what I'm supposed to do and where to find some of these things in the program so I can finish this long and horrible assignment finally.

Here's the problem:

Radionuclides used in PET scanning are typically isotopes with short half-lives such as carbon-11(~20 min), nitrogen-13(~10 min), oxygen-15(~2 min), and fluorine-18(~110 min). You are given 3 mystery materials, which contain different amounts of these four isotopes. You have a dumb counter, which can only measure total counts but cannot discriminate. You measure every minute for 2 hours(121 measurements). Your goal is to determine the amount of each isotope in the 3 mystery containers.

i. Set up least squares problem that you are going to solve. Compute the governing matrix
ii. Plot the columns of the governing matrix vs. time. You should have 4 decaying functions on one plot.
iii. What is the condition number of this matrix?
iv. Solve for each of the noise-free("TRU") datasets
v.Solve for each of the noisy datasets
vi. For which noisy dataset is your result the worst?(Hint: think about singular values involved in each of the fits.) How do they compare?

There's a file she e-mailed to us and it has 6 different spreadsheets for it. I've figured out how to plot it all but I really have no idea what she's asking in the first place or how to use MATLAB to get through it. I've been messing around with this for 4 hours and I've gotten nowhere.

Thanks in advance, here's the link to the .mat file

https://docs.google.com/open?id=0B3KGwKGOPK_OZ29QZ3NsU1JudTQ

Answer 1

[Star Strider]

I suggest you start by considering this a sum of exponentials. You have the isotope half-lives, so you can easily compute the decay constants [ K(n) here ] from them, for example as an anonymous function (with the K(n) vector previously defined in your code so the function can access it from the workspace):

f = @(B,t) B(1).*exp(K(1).*t) + ... + B(4).*exp(K(4).*t);

This can be considered a linear problem. If you also have to estimate the decay constants, then this becomes a nonlinear problem:

f = @(B,t) B(1).*exp(B(2).*t) + ... + B(7).*exp(B(8).*t);

In either situation, with a 121-element data vector, estimating four or eight parameters should not be difficult. You can then compare the decay constants with those you already calculated from the half-lives to see what initial concentrations correspond to what isotopes.

The term ‘governing matrix’ is new to me. I can't even find it on a search. (The best I can guess is that it's the Jacobian, or the variance-covariance matrix of the parameters that can be derived from the Jacobian and the residuals.)

‘CAS’ is also new to me. It doesn't seem to be a sufficiently unique acronym that a search suggests anything specific.

[Walter Roberson]

'CAS' can mean "Computer Algebra System", such as Maple or Mathematica.

Answer 2

[Glenn Niles]

That's one of the problems I'm having. I don't know what a governing matrix is, I e-mailed her over the weekend asking that and what she wants us to do with this problem and she responded with "have you found the governing matrix" skipping over the part of me asking what that even means, since that term is in none of my books or notes.

I'll probably just draw a picture of a moose or something and call it quits since this problem has nothing to do with what we're talking about in class and this is on top of the other 71 problems she assigned that pretty much takes up all my time as is. Frustrations unbound.

Thanks for the help though

[Star Strider]

It's my pleasure to help!

You might be able to get partial credit by performing the regression and estimating the initial concentration coefficients. Here is an example (without added Gaussian noise to the y vector) to get you started:

t = [0:20]';
Iso1 = exp(-0.5*t);       % Decay coefficient calculated from T_1/2 for isotope 1
Iso2 = exp(-2.0*t);       % Decay coefficient calculated from T_1/2 for isotope 2

y = (5*Iso1 + 7*Iso2);

B = [ones(size(Iso1)) Iso1 Iso2]\y

Be sure to set up your data as column vectors, as in the example. (See the documentation for mldivide, and if you have the Statistics Toolbox regress for details.) For noisy data, it is necessary to augment the [Iso1 Iso2] matrix as [ones(size(Iso1)) Iso1 Iso2] in order to get an appropriate fit. (The ones vector generates a parameter analogous to the y-intercept (b) in the typical linear model y=m*x+b.)

I get the impression that ‘governing matrix’ is the Jacobian, very easy to calculate in this situation. It becomes the row vector of the partial derivatives of the function you're fitting with respect to each parameter. It becomes a [121 x 4] matrix when you evaluate it at each time (minute). (If she hasn't defined it in lecture, handouts, or textbook, she can't very well tell you the Jacobian is the wrong answer!) Since she says '4 decaying functions', that implies my first anonymous function (in B and K) describes your problem, however since it is linear (as in the example), you won't actually need to use the function.

Alternatively, the ‘governing matrix’ might be something like:

GM = [Iso1  Iso2] * B(2:end);

I have no idea, but whatever gives you the result closest to what your instructor seems to want is likely the correct one.

Good luck in your course, since it sounds like your instructor is going to be playing ‘guess what I'm thinking’ all term!