| |
 |
|
|
Science Forum Index » Math - Symbolic Forum » Is this a possible Maxima bug or user error?...
Page 1 of 2 Goto page 1, 2 Next
|
| Author |
Message |
| pnachtwey... |
 Posted: Sat Jun 21, 2008 5:54 pm |
|
|
|
Guest
|
Is this a bug or am I doing something wrong?
I am trying to symbolically do a matrix exp function or expm( ). My
Ac is the transition matrix for an underdamped second order system
where the velocity is integrated to a position.
w is the natural frequency
z is the damping factor
Ac: matrix([0,1,0],[0,0,1],[0,-w^2,-2*z*w]);
I: ident(3);
A: I+(Ac*T)+(Ac*T)^2/2!;
grind(A);
The output is:
matrix([1,T^2/2+T,0],[0,1,T^2/2+T], [0,w^4*T^2/2-w^2*T,
2*w^2*z^2*T^2-2*w*z*T+1])$
According to my Mathcad the A[0,2] should ot be zero and A[1,1] should
not be 1.
If I am doing something wrong let me know. I fhis is a bug then it is
very basic.
Peter Nachtwey |
|
|
| Back to top |
|
| rjf... |
| Posted: Sat Jun 21, 2008 7:06 pm |
|
|
|
Guest
|
On Jun 21, 8:54 pm, pnachtwey <pnacht... at (no spam) gmail.com> wrote:
Quote: Is this a bug or am I doing something wrong?
]
You are doing something wrong. if S is a matrix, S^2 squares each
entry. You want S^^2. |
|
|
| Back to top |
|
| pnachtwey... |
| Posted: Sun Jun 22, 2008 7:19 am |
|
|
|
Guest
|
On Jun 21, 10:06 pm, rjf <fate... at (no spam) gmail.com> wrote:
Quote: On Jun 21, 8:54 pm, pnachtwey <pnacht... at (no spam) gmail.com> wrote:> Is this a bug or am I doing something wrong?
]
You are doing something wrong. if S is a matrix, S^2 squares each
entry. You want S^^2.
S.S seems to work too.
S: ((((((M/7+I).M/6+I).M/5+I).M/4+I).M/3+I).M/2+I).M+I;
This works except it is missing the scale and square, after the scale
and square the results will probably be too messy.
I tried doing a
%e^S
%e^^S
exp(S)
just to check my results and the results were all wrong too. I used
Scilab to check my results. It looks like there is a meed for a expm
function like Matlab and Scilab have. I looked all over and didn't
find one. I didn't find a matrix norm either. This is necessary for
the scale and square part of the expm(Matrix). Hmmm. I had to write
this from scratch in Mathcad too.
Thanks
Peter Nachtwey |
|
|
| Back to top |
|
| rjf... |
| Posted: Sun Jun 22, 2008 9:09 am |
|
|
|
Guest
|
On Jun 22, 10:19 am, pnachtwey <pnacht... at (no spam) gmail.com> wrote:
Quote: On Jun 21, 10:06 pm, rjf <fate... at (no spam) gmail.com> wrote:> On Jun 21, 8:54 pm, pnachtwey <pnacht... at (no spam) gmail.com> wrote:> Is this a bug or am I doing something wrong?
]
You are doing something wrong. if S is a matrix, S^2 squares each
entry. You want S^^2.
S.S seems to work too.
yes, that is the same as S^^2.
Quote:
S: ((((((M/7+I).M/6+I).M/5+I).M/4+I).M/3+I).M/2+I).M+I;
This works except it is missing the scale and square, after the scale
and square the results will probably be too messy.
could be why it is not so interesting to do.
Quote:
I tried doing a
%e^S
%e^^S
exp(S)
just to check my results and the results were all wrong too.
I expect that is because Maxima probably doesn't know about this
concept of matrix exponential.
There are many ways of computing a matrix exponential of a numerical
matrix. Except for the approximation by taylor series, I do not know
how to do this in general for symbolic matrices.
(though see
http://www.mathworks.com/access/helpdesk/help/toolbox/symbolic/index.html?/access/helpdesk/help/toolbox/symbolic/expm.html&http://www.google.com/search?q=matrix+exponential+symbolic&sourceid=navclient-ff&ie=UTF-8&rls=GGIC,GGIC:1970--2,GGIC:en
for a special case)
If you have a suggestion for how to do this, |
|
|
| Back to top |
|
| pnachtwey... |
| Posted: Sun Jun 22, 2008 11:09 am |
|
|
|
Guest
|
On Jun 22, 12:22 pm, Axel Vogt <&nore... at (no spam) axelvogt.de> wrote:
Quote: rjf wrote:
On Jun 22, 10:19 am, pnachtwey <pnacht... at (no spam) gmail.com> wrote:
On Jun 21, 10:06 pm, rjf <fate... at (no spam) gmail.com> wrote:> On Jun 21, 8:54 pm, pnachtwey <pnacht... at (no spam) gmail.com> wrote:> Is this a bug or am I doing something wrong?
]
You are doing something wrong. if S is a matrix, S^2 squares each
entry. You want S^^2.
S.S seems to work too.
yes, that is the same as S^^2.
S: ((((((M/7+I).M/6+I).M/5+I).M/4+I).M/3+I).M/2+I).M+I;
This works except it is missing the scale and square, after the scale
and square the results will probably be too messy.
could be why it is not so interesting to do.
I tried doing a
%e^S
%e^^S
exp(S)
just to check my results and the results were all wrong too.
I expect that is because Maxima probably doesn't know about this
concept of matrix exponential.
There are many ways of computing a matrix exponential of a numerical
matrix. Except for the approximation by taylor series, I do not know
how to do this in general for symbolic matrices.
(though see
http://www.mathworks.com/access/helpdesk/help/toolbox/symbolic/index....
for a special case)
If you have a suggestion for how to do this,
There is an overview by Moler & van Loan, "Nineteen Dubious Ways to
Compute the Exponential of a Matrix, Twenty-Five Years Later" (2003)
and their former "nineteen dubious ways" 1978, both SIAM Reviews.- Hide quoted text -
- Show quoted text -
I have read it that document. The Pade approximant with the scale and
square works best but it is very computationally intensive. The
division is a killer. I needed this for a running a simulator for a
type 1 under damped second order system. I am using a DSP and need
6-7 digit of precision.
Ac: matrix([0,1,0],[0,0,1],[0,-w^2,-2*z*w]);
I: ident(3);
A: I+(Ac*T)+(Ac*T)^^2/2!+(Ac*T)^^3/3!+(Ac*T)^^4/4!+(Ac*T)^^5/5!+
(Ac*T)^^6/6!+(Ac*T)^^7/7!;
grind(A);
Ac is the continuous transition matrix
T is the simulation update period.
z is the damping factor
w is the natural frequency.
As you can see Ac is mostly zeros so I was hoping I could reduce the
number of mulitplies and adds by simplifying the math symbolically.
Actually expm(Ac*T) is easy to calculate using Laplace transforms for
2x2 and 3x3 matrices. It is the most accurate and fastest because
only one call to exp, sin and cos is required. Also, there is no need
for the scale and square. Calculating expm() is just a warm
excercise.
The real problem is calculating the discrete time matrix for B in
x(n+1)=A*x(n)+B*u(n);
This equation is used in simulation.
Bc: matrix( [0],[0],[K*w^2]); /* Input coupling matrix */
K is the system gain
B: (I+(Ac*T)/2!+(Ac*T)^^2/3!+(Ac*T)^^3/4!+(Ac*T)^^4/5!+(Ac*T)^^5/6!+
(Ac*T)^^6/7!).(Bc*T);
Compare the symbolic result with calculating with Taylor's method even
when calculated more efficiently.
S: ((((((M/7+I).M/6+I).M/5+I).M/4+I).M/3+I).M/2+I).M+I;
Since only three values are required why not calculate then directly?
This still leaves the scale and square issue. If the damping factor
and natural frquency are at the higher end of the range the values in
the lower right part of the A matrix will not converge quickly. The
scale and square technique is required.
ftp://ftp.deltacompsys.com/public/NG/Mathcad%20-%20expm.pdf
You can see I have done this before. After page 1 it is just random
calculations. I posted this some time back to show how to calcuate
the expm() but I was using Mathcad then. Mathcad does not let me cut
and paste equations into C. Maxima does but I most edit a little.
Peter Nachtwey |
|
|
| Back to top |
|
| Jim... |
| Posted: Sun Jun 22, 2008 11:51 am |
|
|
|
Guest
|
Sun, 22 Jun 2008 21:22:05 +0200, Axel Vogt:
Quote: rjf wrote:
I expect that is because Maxima probably doesn't know about this
concept of matrix exponential.
Nonsense.
A : matrix([2,1],[0,2]);
load(diag);
mat_function(exp,A);
matrix([%e^2,%e^2],[0,%e^2]) |
|
|
| Back to top |
|
| rjf... |
| Posted: Sun Jun 22, 2008 1:48 pm |
|
|
|
Guest
|
On Jun 22, 2:51 pm, Jim <J... at (no spam) j.eu> wrote:
Quote: Sun, 22 Jun 2008 21:22:05 +0200, Axel Vogt:
rjf wrote:
I expect that is because Maxima probably doesn't know about this
concept of matrix exponential.
Nonsense.
A : matrix([2,1],[0,2]);
load(diag);
mat_function(exp,A);
matrix([%e^2,%e^2],[0,%e^2])
Well, I stand corrected. Thanks. Maxima knows something about matrix
exponentials (and other analytic functions,... you can do
mat_function(cos,A) also).
On the other hand, beyond 2X2 or even general 2X2, it is not going to
do anything magic.
You might be lucky, but e.g. if v is a general 5X5 symbolic matrix,
mat_function(exp,v) says "solve is unable to find the roots of the
characteristic polynomial...."
This is a mathematical difficulty, not a defect of Maxima
particularly... |
|
|
| Back to top |
|
| Axel Vogt... |
| Posted: Sun Jun 22, 2008 2:22 pm |
|
|
|
Guest
|
rjf wrote:
Quote: On Jun 22, 10:19 am, pnachtwey <pnacht... at (no spam) gmail.com> wrote:
On Jun 21, 10:06 pm, rjf <fate... at (no spam) gmail.com> wrote:> On Jun 21, 8:54 pm, pnachtwey <pnacht... at (no spam) gmail.com> wrote:> Is this a bug or am I doing something wrong?
]
You are doing something wrong. if S is a matrix, S^2 squares each
entry. You want S^^2.
S.S seems to work too.
yes, that is the same as S^^2.
S: ((((((M/7+I).M/6+I).M/5+I).M/4+I).M/3+I).M/2+I).M+I;
This works except it is missing the scale and square, after the scale
and square the results will probably be too messy.
could be why it is not so interesting to do.
I tried doing a
%e^S
%e^^S
exp(S)
just to check my results and the results were all wrong too.
I expect that is because Maxima probably doesn't know about this
concept of matrix exponential.
There are many ways of computing a matrix exponential of a numerical
matrix. Except for the approximation by taylor series, I do not know
how to do this in general for symbolic matrices.
(though see
http://www.mathworks.com/access/helpdesk/help/toolbox/symbolic/index.html?/access/helpdesk/help/toolbox/symbolic/expm.html&http://www.google.com/search?q=matrix+exponential+symbolic&sourceid=navclient-ff&ie=UTF-8&rls=GGIC,GGIC:1970--2,GGIC:en
for a special case)
If you have a suggestion for how to do this,
There is an overview by Moler & van Loan, "Nineteen Dubious Ways to
Compute the Exponential of a Matrix, Twenty-Five Years Later" (2003)
and their former "nineteen dubious ways" 1978, both SIAM Reviews. |
|
|
| Back to top |
|
| Jim... |
| Posted: Mon Jun 23, 2008 2:59 am |
|
|
|
Guest
|
Sun, 22 Jun 2008 16:48:46 -0700, rjf:
Quote: On Jun 22, 2:51 pm, Jim <J... at (no spam) j.eu> wrote:
Sun, 22 Jun 2008 21:22:05 +0200, Axel Vogt:
rjf wrote:
I expect that is because Maxima probably doesn't know about this
concept of matrix exponential.
Nonsense.
A : matrix([2,1],[0,2]);
load(diag);
mat_function(exp,A);
matrix([%e^2,%e^2],[0,%e^2])
Well, I stand corrected. Thanks. Maxima knows something about matrix
exponentials (and other analytic functions,... you can do
mat_function(cos,A) also).
On the other hand, beyond 2X2 or even general 2X2, it is not going to do
anything magic.
You might be lucky, but e.g. if v is a general 5X5 symbolic matrix,
mat_function(exp,v) says "solve is unable to find the roots of the
characteristic polynomial...."
How strange... Galois should know something about it...
Jokes apart, on a generic symbolic matrix it's fairly useless:
M: matrix([a,b],[c,d])
mat_function(exp,M);
$$\pmatrix{\frac{\left( \sqrt{{d}^{2}-2\,a\,d+4\,b\,c+{a}^{2}}-d+a\right)
\,{e}^{\frac{\sqrt{{d}^{2}-2\,a\,d+4\,b\,c+{a}^{2}}+d+a}{2}}}{2\,\sqrt{{d}^
{2}-2\,a\,d+4\,b\,c+{a}^{2}}}+\frac{\left( \sqrt{{d}^{2}-2\,a\,d+4\,b\,c
+{a}^{2}}+d-a\right) \,{e}^{-\frac{\sqrt{{d}^{2}-2\,a\,d+4\,b\,c+{a}^{2}}-
d-a}{2}}}{2\,\sqrt{{d}^{2}-2\,a\,d+4\,b\,c+{a}^{2}}} & \frac{b\,{e}^{\frac
{\sqrt{{d}^{2}-2\,a\,d+4\,b\,c+{a}^{2}}+d+a}{2}}}{\sqrt{{d}^{2}-2\,a\,d+4
\,b\,c+{a}^{2}}}-\frac{b\,{e}^{-\frac{\sqrt{{d}^{2}-2\,a\,d+4\,b\,c+{a}^
{2}}-d-a}{2}}}{\sqrt{{d}^{2}-2\,a\,d+4\,b\,c+{a}^{2}}}\cr \frac{c\,{e}^
{\frac{\sqrt{{d}^{2}-2\,a\,d+4\,b\,c+{a}^{2}}+d+a}{2}}}{\sqrt{{d}^{2}-2\,a
\,d+4\,b\,c+{a}^{2}}}-\frac{c\,{e}^{-\frac{\sqrt{{d}^{2}-2\,a\,d+4\,b\,c
+{a}^{2}}-d-a}{2}}}{\sqrt{{d}^{2}-2\,a\,d+4\,b\,c+{a}^{2}}} & \frac{2\,b\,c
\,{e}^{\frac{\sqrt{{d}^{2}-2\,a\,d+4\,b\,c+{a}^{2}}+d+a}{2}}}{\sqrt{{d}^
{2}-2\,a\,d+4\,b\,c+{a}^{2}}\,\left( \sqrt{{d}^{2}-2\,a\,d+4\,b\,c+{a}^
{2}}-d+a\right) }+\frac{2\,b\,c\,{e}^{-\frac{\sqrt{{d}^{2}-2\,a\,d+4\,b\,c
+{a}^{2}}-d-a}{2}}}{\sqrt{{d}^{2}-2\,a\,d+4\,b\,c+{a}^{2}}\,\left( \sqrt
{{d}^{2}-2\,a\,d+4\,b\,c+{a}^{2}}+d-a\right) }}$$ |
|
|
| Back to top |
|
| rjf... |
| Posted: Mon Jun 23, 2008 5:11 am |
|
|
|
Guest
|
On Jun 23, 5:59 am, Jim <J... at (no spam) j.eu> wrote:
Quote: load(diag);
mat_function(exp,A);
matrix([%e^2,%e^2],[0,%e^2])
Well, I stand corrected. Thanks. Maxima knows something about matrix
exponentials (and other analytic functions,... you can do
mat_function(cos,A) also).
On the other hand, beyond 2X2 or even general 2X2, it is not going to do
anything magic.
You might be lucky, but e.g. if v is a general 5X5 symbolic matrix,
mat_function(exp,v) says "solve is unable to find the roots of the
characteristic polynomial...."
How strange... Galois should know something about it...
Jokes apart, on a generic symbolic matrix it's fairly useless:
M: matrix([a,b],[c,d])...
<snip>
on the other hand, you can try the Maxima program for matrix
exponential on a diagonal matrix of any size and get very satisfactory
results.
RJF |
|
|
| Back to top |
|
| Jim... |
| Posted: Mon Jun 23, 2008 6:11 am |
|
|
|
Guest
|
Mon, 23 Jun 2008 08:11:23 -0700, rjf:
Quote: on the other hand, you can try the Maxima program for matrix exponential
on a diagonal matrix of any size and get very satisfactory results.
RJF
Why should one do that? |
|
|
| Back to top |
|
| pnachtwey... |
| Posted: Mon Jun 23, 2008 11:56 am |
|
|
|
Guest
|
On Jun 22, 4:48 pm, rjf <fate... at (no spam) gmail.com> wrote:
Quote: On Jun 22, 2:51 pm, Jim <J... at (no spam) j.eu> wrote:
Sun, 22 Jun 2008 21:22:05 +0200, Axel Vogt:
rjf wrote:
I expect that is because Maxima probably doesn't know about this
concept of matrix exponential.
Nonsense.
A : matrix([2,1],[0,2]);
load(diag);
mat_function(exp,A);
matrix([%e^2,%e^2],[0,%e^2])
Well, I stand corrected. Thanks. Maxima knows something about matrix
exponentials (and other analytic functions,... you can do
mat_function(cos,A) also).
On the other hand, beyond 2X2 or even general 2X2, it is not going to
do anything magic.
You might be lucky, but e.g. if v is a general 5X5 symbolic matrix,
mat_function(exp,v) says "solve is unable to find the roots of the
characteristic polynomial...."
This is a mathematical difficulty, not a defect of Maxima
particularly...
I didn't have much luck with the mat_function(exp,mat). I couldn't
get it to work except on the simplest of arrays as shown above.
Since you guys are being so helpful I have a different approach to the
expm(Ac*T) problem.
Ac: matrix([0,1,0],[0,0,1],[0,-w^2,-2*z*w]);
I: ident(3);
q0: (s*I-Ac)^^-1;
A[1,1]: ilt(q0[1,1],s,t); /* since z is < 1 the answer is no */
A[1,2]: ilt(q0[1,2],s,t);
A[1,3]: ilt(q0[1,3],s,t); /* This one will not work because of the
questions */
A[2,1]: ilt(q0[2,1],s,t);
A[2,2]: ilt(q0[2,2],s,t);
A[2,3]: ilt(q0[2,3],s,t);
A[3,1]: ilt(q0[3,1],s,t);
A[3,2]: ilt(q0[3,2],s,t);
A[3,3]: ilt(q0[3,3],s,t);
I can do this in Mathcad by doing a lot of cutting and pasting. See
the link I posted to the expm.pdf. I calculated expm(Ac*T) using
Laplaces transforms. I was wondering if Maxima had a better way where
these calculations can be calculated with little manual interaction.
It would be nice if I could just:
A: ilt(q0,s,t)
There seem to be a couple of problems. First, Maxima asks me a
question for which I must type an answer. This does not help in
automating things but it is no problem if I want to answer the same
question nxn times. It would be nice to store the invlaplace in an
array A so all I need to do is provide values for T,z and w and out
pops a discrete transition matrix expm(Ac*T). I tried using nested
for loops but I think the nested for loops dont' work because of
question that is repeatedly asksed. Is there a way of have a default
answer for the question since it is the same each time?
This is my prefered method to calculate expm since I don't need to
worry about the scale and square and the calculation are really quite
fast because there are many common terms between the different cells
in A.
Can Maple or Mathematica do these kinds of problems easier? I keep
getting burned by Maxima's odd syntax, poor documentation and lack of
examples.
Hopefully some of you are getting something out of this.
Thanks
Peter Nachtwey |
|
|
| Back to top |
|
| Jim... |
| Posted: Mon Jun 23, 2008 1:02 pm |
|
|
|
Guest
|
Mon, 23 Jun 2008 14:56:21 -0700, pnachtwey:
Quote: I didn't have much luck with the mat_function(exp,mat). I couldn't get
it to work except on the simplest of arrays as shown above.
It works for your Ac matrix, see below.
Since you guys are being so helpful I have a different approach to the
expm(Ac*T) problem.
Ac: matrix([0,1,0],[0,0,1],[0,-w^2,-2*z*w]);
mat_function(exp,Ac);
matrix([1,((2*sqrt(z-1)*z*sqrt(z+1)+2*z^2-1)*%e^(w*sqrt(z^2-1)-w*z))/(sqrt
(z^2-1)*(4*w*z^2-2*w)-4*w*sqrt(z-1)*z^2*sqrt(z+1))+((2*sqrt(z-1)*z*sqrt(z
+1)-2*z^2+1)*%e^(-w*sqrt(z^2-1)-w*z))/(sqrt(z^2-1)*(4*w*z^2-2*w)-4*w*sqrt
(z-1)*z^2*sqrt(z+1))-(2*sqrt(z-1)*z*sqrt(z+1))/(sqrt(z^2-1)*(2*w*z^2-
w)-2*w*sqrt(z-1)*z^2*sqrt(z+1)),((sqrt(z^2-1)+z)*%e^(w*sqrt(z^2-1)-w*z))/
(sqrt(z^2-1)*(4*w^2*z^2-2*w^2)-4*w^2*sqrt(z-1)*z^2*sqrt(z+1))+((sqrt
(z^2-1)-z)*%e^(-w*sqrt(z^2-1)-w*z))/(sqrt(z^2-1)*
(4*w^2*z^2-2*w^2)-4*w^2*sqrt(z-1)*z^2*sqrt(z+1))-sqrt(z^2-1)/(sqrt(z^2-1)*
(2*w^2*z^2-w^2)-2*w^2*sqrt(z-1)*z^2*sqrt(z+1))],[0,((2*sqrt(z-1)*z*sqrt(z
+1)+2*z^2-1)*(w*sqrt(z^2-1)-w*z)*%e^(w*sqrt(z^2-1)-w*z))/(sqrt(z^2-1)*
(4*w*z^2-2*w)-4*w*sqrt(z-1)*z^2*sqrt(z+1))+((2*sqrt(z-1)*z*sqrt(z+1)-2*z^2
+1)*(-w*sqrt(z^2-1)-w*z)*%e^(-w*sqrt(z^2-1)-w*z))/(sqrt(z^2-1)*
(4*w*z^2-2*w)-4*w*sqrt(z-1)*z^2*sqrt(z+1)),((sqrt(z^2-1)+z)*(w*sqrt(z^2-1)-
w*z)*%e^(w*sqrt(z^2-1)-w*z))/(sqrt(z^2-1)*(4*w^2*z^2-2*w^2)-4*w^2*sqrt(z-1)
*z^2*sqrt(z+1))+((sqrt(z^2-1)-z)*(-w*sqrt(z^2-1)-w*z)*%e^(-w*sqrt(z^2-1)-
w*z))/(sqrt(z^2-1)*(4*w^2*z^2-2*w^2)-4*w^2*sqrt(z-1)*z^2*sqrt(z+1))],[0,
((2*sqrt(z-1)*z*sqrt(z+1)+2*z^2-1)*(-2*w^2*sqrt(z-1)*z*sqrt(z+1)+2*w^2*z^2-
w^2)*%e^(w*sqrt(z^2-1)-w*z))/(sqrt(z^2-1)*(4*w*z^2-2*w)-4*w*sqrt(z-1)
*z^2*sqrt(z+1))+((2*sqrt(z-1)*z*sqrt(z+1)-2*z^2+1)*(2*w^2*sqrt(z-1)*z*sqrt
(z+1)+2*w^2*z^2-w^2)*%e^(-w*sqrt(z^2-1)-w*z))/(sqrt(z^2-1)*
(4*w*z^2-2*w)-4*w*sqrt(z-1)*z^2*sqrt(z+1)),((-2*w^2*sqrt(z-1)*z*sqrt(z+1)
+2*w^2*z^2-w^2)*(sqrt(z^2-1)+z)*%e^(w*sqrt(z^2-1)-w*z))/(sqrt(z^2-1)*
(4*w^2*z^2-2*w^2)-4*w^2*sqrt(z-1)*z^2*sqrt(z+1))+((2*w^2*sqrt(z-1)*z*sqrt(z
+1)+2*w^2*z^2-w^2)*(sqrt(z^2-1)-z)*%e^(-w*sqrt(z^2-1)-w*z))/(sqrt(z^2-1)*
(4*w^2*z^2-2*w^2)-4*w^2*sqrt(z-1)*z^2*sqrt(z+1))]) |
|
|
| Back to top |
|
| Jim... |
| Posted: Mon Jun 23, 2008 1:05 pm |
|
|
|
Guest
|
Mon, 23 Jun 2008 23:02:21 +0000, Jim:
Quote: Ac: matrix([0,1,0],[0,0,1],[0,-w^2,-2*z*w]);
And you need to plug your time evolution, right?
mat_function(exp,Ac*t);
matrix([1,((2*sqrt(z-1)*z*sqrt(z+1)+2*z^2-1)*%e^(t*w*sqrt(z^2-1)-t*w*z))/
(sqrt(z^2-1)*(4*w*z^2-2*w)-4*w*sqrt(z-1)*z^2*sqrt(z+1))+((2*sqrt(z-1)
*z*sqrt(z+1)-2*z^2+1)*%e^(-t*w*sqrt(z^2-1)-t*w*z))/(sqrt(z^2-1)*
(4*w*z^2-2*w)-4*w*sqrt(z-1)*z^2*sqrt(z+1))-(2*sqrt(z-1)*z*sqrt(z+1))/(sqrt
(z^2-1)*(2*w*z^2-w)-2*w*sqrt(z-1)*z^2*sqrt(z+1)),((sqrt(z^2-1)+z)*%e^
(t*w*sqrt(z^2-1)-t*w*z))/(sqrt(z^2-1)*(4*w^2*z^2-2*w^2)-4*w^2*sqrt(z-1)
*z^2*sqrt(z+1))+((sqrt(z^2-1)-z)*%e^(-t*w*sqrt(z^2-1)-t*w*z))/(sqrt(z^2-1)*
(4*w^2*z^2-2*w^2)-4*w^2*sqrt(z-1)*z^2*sqrt(z+1))-sqrt(z^2-1)/(sqrt(z^2-1)*
(2*w^2*z^2-w^2)-2*w^2*sqrt(z-1)*z^2*sqrt(z+1))],[0,((2*sqrt(z-1)*z*sqrt(z
+1)+2*z^2-1)*(w*sqrt(z^2-1)-w*z)*%e^(t*w*sqrt(z^2-1)-t*w*z))/(sqrt(z^2-1)*
(4*w*z^2-2*w)-4*w*sqrt(z-1)*z^2*sqrt(z+1))+((2*sqrt(z-1)*z*sqrt(z+1)-2*z^2
+1)*(-w*sqrt(z^2-1)-w*z)*%e^(-t*w*sqrt(z^2-1)-t*w*z))/(sqrt(z^2-1)*
(4*w*z^2-2*w)-4*w*sqrt(z-1)*z^2*sqrt(z+1)),((sqrt(z^2-1)+z)*(w*sqrt(z^2-1)-
w*z)*%e^(t*w*sqrt(z^2-1)-t*w*z))/(sqrt(z^2-1)*(4*w^2*z^2-2*w^2)-4*w^2*sqrt
(z-1)*z^2*sqrt(z+1))+((sqrt(z^2-1)-z)*(-w*sqrt(z^2-1)-w*z)*%e^(-t*w*sqrt
(z^2-1)-t*w*z))/(sqrt(z^2-1)*(4*w^2*z^2-2*w^2)-4*w^2*sqrt(z-1)*z^2*sqrt(z
+1))],[0,((2*sqrt(z-1)*z*sqrt(z+1)+2*z^2-1)*(-2*w^2*sqrt(z-1)*z*sqrt(z+1)
+2*w^2*z^2-w^2)*%e^(t*w*sqrt(z^2-1)-t*w*z))/(sqrt(z^2-1)*
(4*w*z^2-2*w)-4*w*sqrt(z-1)*z^2*sqrt(z+1))+((2*sqrt(z-1)*z*sqrt(z+1)-2*z^2
+1)*(2*w^2*sqrt(z-1)*z*sqrt(z+1)+2*w^2*z^2-w^2)*%e^(-t*w*sqrt(z^2-1)-
t*w*z))/(sqrt(z^2-1)*(4*w*z^2-2*w)-4*w*sqrt(z-1)*z^2*sqrt(z+1)),
((-2*w^2*sqrt(z-1)*z*sqrt(z+1)+2*w^2*z^2-w^2)*(sqrt(z^2-1)+z)*%e^(t*w*sqrt
(z^2-1)-t*w*z))/(sqrt(z^2-1)*(4*w^2*z^2-2*w^2)-4*w^2*sqrt(z-1)*z^2*sqrt(z
+1))+((2*w^2*sqrt(z-1)*z*sqrt(z+1)+2*w^2*z^2-w^2)*(sqrt(z^2-1)-z)*%e^(-
t*w*sqrt(z^2-1)-t*w*z))/(sqrt(z^2-1)*(4*w^2*z^2-2*w^2)-4*w^2*sqrt(z-1)
*z^2*sqrt(z+1))]) |
|
|
| Back to top |
|
| Christopher Creutzig... |
| Posted: Mon Jun 23, 2008 2:35 pm |
|
|
|
Guest
|
rjf wrote:
Quote: There are many ways of computing a matrix exponential of a numerical
matrix. Except for the approximation by taylor series, I do not know
how to do this in general for symbolic matrices.
Certainly not in general, but for the handful of symbolic matrices with
a computable Jordan form, you can compute the exact value of the Taylor
series.
--
Quote: seit wann sind Vertragsinhalte für NewsGroup-Frager relevant?
Sie sind lebensnotwendig um sie sofort auf überraschende Inhalte
abzuklopfen oder sonst in Frage zu stellen.
(Kurt Gunter und Konrad Wilhelm in dsrm) |
|
|
| Back to top |
|
| |
Page 1 of 2 Goto page 1, 2 Next
All times are GMT - 5 Hours
The time now is Fri Dec 05, 2008 9:46 am
|
|