Message Boards Message Boards

0
|
6753 Views
|
11 Replies
|
2 Total Likes
View groups...
Share
Share this post:

Distributed parameter beam system: determine natural frequencies

Kaitlin Smith
Posted 11 years ago
Hi,

I am working on a method to evaluate a distributed parameter beam system to determine its natural frequencies.  To do so, I need to find the exponential matrix of a six by six symbolic matrix, and then integrate the individual entries using the symbol as the variable from 0 to 1.

So I start with this matrix:
bigPhis = {{0, 1, 0, 0, 0, 0}, {0, 0, 1, 0, 0, 0}, {0, 0, 0, 1, 0,
    0}, {0, 0, 0, 0, 1, 0}, {0, 0, 0, 0, 0, 1}, {c1, 0, c2, 0,
    alpha^2, 0}};

where c1 is a complex number, and c2 and alpha are both real.  I am able to get the exponential matrix output (point for Mathematica over Matlab!) using
bigPsixs = MatrixExp[bigPhis*(x - xo)]

where x is a symbol and xo is a value of 0.5.  I take the first entry of this matrix, which gives me an output that starts like this:


{RootSum[(8.29521*10^6 -
      55.5556 I) - (9.95425*10^7 - 666.667 I) x + (4.97712*10^8 -
       3333.33 I) x^2 - (1.32723*10^9 -
       8888.89 I) x^3 + (1.99085*10^9 -
       13333.3 I) x^4 - (1.59268*10^9 -
       10666.7 I) x^5 + (5.30893*10^8 - 3555.56 I) x^6 -
    331809. #1^2 + 3.98171*10^6 x #1^2 - 1.99086*10^7 x^2 #1^2 +
    5.30895*10^7 x^3 #1^2 - 7.96342*10^7 x^4 #1^2 + ... ...  ...
(continues)


So I have read the articles about #1 being a placeholder for the first root solution, which I kind of sort of understand, but the next step for this problem is:

om11 = Integrate[Exp[alpha*(x - xo)]*bigPsixs[[1, 1]], {x, 0, xo}] +  Integrate[Exp[alpha*(xo - x)]*bigPsixs[[1, 1]], {x, xo, 1}];

So you can see that I still need x in the bigPsixs input to be able to integrate it.  Does anyone have any tips or ideas of how to make this work?  I don't really know how to solve for #1 numerically and still keep x in the equation for the next step integration, but maybe I am confused on the Mathematica output... please help?

Thanks,
Kaitlin
POSTED BY: Kaitlin Smith
11 Replies
Sort By:
Probably better to use NIntegrate for this. Alternatively, use exact values everywhere and do
Integrate[Exp[1/1000*(x - 1/2)]*bigPsixs[[1, 1]], {x, 0, 5/100}];
It seems much faster to avoid the ToRadicals[...] stuff.
POSTED BY: Daniel Lichtblau
Bill Simpson
Bill Simpson
Posted 11 years ago
I assume you are using the notebook interface. I assume you executed the example I showed with your expression pasted into the RootSum. I assume you then see a box with "A very large output was generated. Here is a sample of it:" followed by 0.16667 (<<1>>) followed by a button that says "Show Full Output." and if you click on that then you should see an expression that begins with 0.16667(( followed by about 200 screenfuls that is all the ugly details of the expression you just asked for. That result appears to be purely in terms of x and Euler's constant and powers and radicals and all the usual algebraic and transcendental things you expect in a solution, the RootSum and #1 have all been removed in the half dozen screens that I looked at. Searching the output for a # wraps around and only finds those in your original expression, so I believe this is strong evidence that this did as you asked, got rid of the RootSum and #1. If you wrap an N[] around the ToRadicals then it does decimal approximations everwhere it can and "Show Full Output" shows this reduced down to about 100 screenfulls. Verify that you see both those.

If you then expect to integrate that a few times then my best guess would be that this will not finish in your lifetime.

Can you go back the the underlying pure mathematics behind your problem and find a way to make the problem hundreds or thousands of times simpler?

Or use Daniel's brillian suggestion that he offered while I was typing all this.
POSTED BY: Bill Simpson
The symbol name fell out of "the symbol as the variable".
POSTED BY: Bruce Miller
Kaitlin Smith
Kaitlin Smith
Posted 11 years ago
Sorry, I am not quite sure what you mean- "the symbol as the variable" is x- is that what you are referring to?  Why or how would a variable "fall out" and how can I fix it?
POSTED BY: Kaitlin Smith
Bill Simpson
Bill Simpson
Posted 11 years ago
In[1]:= bigPhis = {{0,1,0,0,0,0}, {0,0,1,0,0,0}, {0,0,0,1,0,
0}, {0,0,0,0,1,0}, {0,0,0,0,0,1}, {c1,0,c2,0,alpha^2,0}};
bigPsixs = MatrixExp[bigPhis*(x - xo)];
ToRadicals[bigPsixs[[1]]]

Out[3]= ...huge output with x snipped, perhaps with your coefficients it will be smaller, perhaps Simplify may help
POSTED BY: Bill Simpson
Kaitlin Smith
Kaitlin Smith
Posted 11 years ago
I like the ToRadicals idea, but I don't think it worked- here is what I got from "testy=ToRadicals[bigPsixs[[1,1]]" (I only need one entry at a time, so I changed [[1]] to [[1,1]]).
I still end up with #1 in the output...
 RootSum[(0.0801746 -
      0.00555556 I) - (0.962095 - 0.0666667 I) x + (4.81047 -
       0.333333 I) x^2 - (12.8279 - 0.888889 I) x^3 + (19.2419 -
       1.33333 I) x^4 - (15.3935 - 1.06667 I) x^5 + (5.13117 -
       0.355556 I) x^6 - 331809. #1^2 + 2.65447*10^6 x #1^2 -
    7.96342*10^6 x^2 #1^2 + 1.06179*10^7 x^3 #1^2 -
    5.30895*10^6 x^4 #1^2 - 2.5*10^-7 #1^4 + 1.*10^-6 x #1^4 -
    1.*10^-6 x^2 #1^4 +
    1. #1^6 &, ((0. + 0. I) - 331809. E^#1 #1 +
     2.65447*10^6 E^#1 x #1 - 7.96342*10^6 E^#1 x^2 #1 +
     1.06179*10^7 E^#1 x^3 #1 - 5.30895*10^6 E^#1 x^4 #1 -
     2.5*10^-7 E^#1 #1^3 + 1.*10^-6 E^#1 x #1^3 -
     1.*10^-6 E^#1 x^2 #1^3 + E^#1 #1^5)/(-663619. #1 +
     5.30895*10^6 x #1 - 1.59268*10^7 x^2 #1 + 2.12358*10^7 x^3 #1 -
     1.06179*10^7 x^4 #1 - 1.*10^-6 #1^3 + 4.*10^-6 x #1^3 -
     4.*10^-6 x^2 #1^3 + 6. #1^5) &]
The coefficients for the code are
c1=-5.13117 + 0.355556 I
c2=5.30895*10^6
alpha=0.001

if you want to play with it using those.  Likely to be no reduction in the output, sorry!  Anyone have ideas on how to eliminate/solve #1 so that I have all of the output in terms of x???
POSTED BY: Kaitlin Smith
Dear Kaitlin, please take a few minutes to read this tutorial about correct posting – especially of Mathematica code:

How to type up a post: editor tutorial & general tips

If you will not follow the tutorial, other members of community may not be able to test your code. To edit your post – click “Edit” in the lower right corner of your post.
POSTED BY: EDITORIAL BOARD
Bill Simpson
Bill Simpson
Posted 11 years ago
c1 = -5.13117 + 0.355556 I;
c2 = 5.30895*10^6;
alpha = 0.001;
ToRadicals[Normal[RootSum[...]]]
with or without an N[] wrapped around the ToRadicals seems to eliminate the #1 and give you the x you desire..
The use of Normal is described in the help page for RootSum hidden under Details
POSTED BY: Bill Simpson
Kaitlin Smith
Kaitlin Smith
Posted 11 years ago
Interesting... that does seem to be an improvement, but is <<1>> that I see in the output still a root solution?  I am trying to get to a point of being able to evaluate this:
om11 = Integrate[
   Exp[0.001*(x - 0.5)]*ToRadicals[Normal[bigPsixs[[1, 1]]]], {x, 0,
    0.05}] + Integrate[
   Exp[0.001*(0.5 - x)]*ToRadicals[Normal[bigPsixs[[1, 1]]]], {x, 0.05,
    1}]

where the bigPsixs matrix is:

bigPhis = {{0, 1, 0, 0, 0, 0}, {0, 0, 1, 0, 0, 0}, {0, 0, 0, 1, 0,
    0}, {0, 0, 0, 0, 1, 0}, {0, 0, 0, 0, 0, 1}, {-5.13117 + 0.355556 I, 0, 5.30895*10^6, 0,
    0.001^2, 0}};
bigPsixs = MatrixExp[bigPhis*(x - 0.5)];

When I try to evaluate even just the first element of bigPsixs by itself, it just runs and runs- whether I use N[ or not.  I had initially assumed it couldn't calculate the integral because I had the #1 terms in the bigPsixs expression, which now seem to be replaced with <<num>> values.  How do I know whether the integration isn't completing because it is computationally intensive or whether it isn't completing because there is a problem with the bigPsixs output?  I am fairly new to Mathamatica, so I don't know how to deal with nasty math very well...

Integrate[ToRadicals[Normal[bigPsixs[[1, 1]]]], {x, 0, 1}]
Integrate[N[ToRadicals[Normal[bigPsixs[[1, 1]]]]], {x, 0, 1}]

Thank you for your help!
POSTED BY: Kaitlin Smith
Kaitlin Smith
Kaitlin Smith
Posted 11 years ago
Daniel Lichtblau wins!  NIntegrate gave me a value within seconds- will now start expanding and see if it can keep up with the rest... thank you all for your suggestions!
POSTED BY: Kaitlin Smith
Kaitlin Smith
Kaitlin Smith
Posted 11 years ago
Reply to Bill Simpson: Yes, those are correct on the outputs- and yes, I totally agree with you that it does not seem to be a great way to approach the problem, I am following a method from a schoalrly paper and am in contact with the authors- unfortunately, the only one who hasn't replied is the one who actually did the problem, and the other two don't seem to remember how it was done- eek!
POSTED BY: Kaitlin Smith
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Tag limit exceeded
Note: Only the first five people you tag will receive an email notification; the other tagged names will appear as links to their profiles.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract