I have come across an integral involving matrices while trying to develop a closed form solution for the problem I have been working on.
I am quite confident about my closed form expression. However, my closed form solution differs from solution obtained using simulation. So, I believe that I am wrong in evaluating that integral.
The integral that I am evaluating at t= 0.0300 sec (end point of Gaussian pulse :
I = Integral (I1. dt)
Where I1 = exp(-I0) ; I0 = Integral((A_mo + Alpha_mo * G(t)) dt)
Both A_mo and Alpha_mo are constant matrices (6 x 6) and G(t) is a Gaussian function.
My Matlab code is attached.
Any help would be highly appreciated.
Do the matrices A_mo and Alpha_mo commute? That is, does A_mo * Alpha_mo = Alpha_mo * A_mo? If not, then the "closed form" solution should be re-checked, as it may be different to a closed form value for the case where A_mo and Alpha_mo are scalars. How did you get the "closed form" solution?
Professor Stewart,
Thanks for responding.
Matrices A-mo and Alpha_mo do NOT commute and I have not commuted them like scalar.
This problem is arising in a MRI application where I am trying to explain some novel contrast. I can easily derive it using simulation. But, then I have to invert my data as well. Since simulation is almost 800 times slower, I needed a closed form solution. However, there are some discrepancies between "simulation" and "closed form solution" and that tells me that there is something wrong in closed form solution.
Based on my investigations, I am reasonably certain that my numerical evaluation of integral is wrong.
Regards,
DK
Since the integral is a linear operator, we have
I0 = Integral((A_mo + Alpha_mo * G(t)) dt)
I0 =A_mo* Integral(dt) + Alpha_mo * Integral(G(t) dt)
If the integral is from 0 to t, say, we have
I0 =A_mo* t + Alpha_mo * Integral(G(t)) dt).
The remaining integral here is essentially an error function.
Thus, you can do the innermost integral analytically. The result is a matrix.
I noticed that A_mo contains elements larger than 1000, leading to elements as large as 30 in A_mo*t for t=3/100 (in absolute value). I suggest that you check whether the expm function can handle such values without underflow or overflow.
Dr. Homeier,
Thanks for responding.
Of course, the evaluation of integral I0 is not a problem.
I believe that somehow my numerical evaluation of: I = Integral (I1. dt) is not correct.
Also, I do not believe that evaluation of expm(…) is the problem as the closed form expression for rectangular pulse (i.e. G(t) = 1 everywhere during the pulse) is matching well with simulation based solution.
Would you be kind enough to look at my derivation and original code as it seems that I am almost out of my depth (I am an experimental physicist!)?
I would be more than happy to invite you to be my co-author upon the publication of my work based on this collaboration.
Regards,
Dushyant
- Badges
- Science topic