How to solve for eigenvalues

Hi,

I am currently investigating a spring restrained, hinged free beam in order to obtain its eigenvalues.

The characteristic equation is: [(αl/k)*[{sin(αl)*cosh(αl)}-{cos(αl)*sinh(αl)}]]=1-(cos(αl)*cosh(αl)).

Where the eigenvalues are αl. A numerical iteration must be done to obtain these eigenvalues, however I am uncertain of how to run this on matlab as I am new to this. Also the value of k can be chosen, so for this I chose 0.5.

I am unsure that I have explained this sufficiently but I can provide anymore information if required. Any help would be greatly appreciated.

Many Thanks,

Bastian Harvey

Answer 1

[John D'Errico]

Learn to use a rootfinder. Also, learn to use valid syntax. {} is not a valid class of parens in a mathematical expression in MATLAB. Also note the use of the .* and ./ operators as I used them here:

eq = @(al,k) ((al/k).*((sin(al).*cosh(al))-(cos(al).*sinh(al)))) - (1-(cos(al).*cosh(al)));
ezplot(@(k) eq(k,0.5),[0,50])
grid on

Your function is symmetric around zero, so I plotted only the positive branch. It appears there are multiple roots, probably infinitely many. If we plot it closer to zero, we can see that zero is a root, as well as a root near 4.

k = 0.5;
syms al
EQ = ((al/k).*((sin(al).*cosh(al))-(cos(al).*sinh(al)))) - (1-(cos(al).*cosh(al)));

The various plots I did indicated a few solutions. You can see some of them here:

vpasolve(EQ,al,4)
ans =
3.8530504836173137215418129889578

vpasolve(EQ,al,20)
ans =
19.622049610416241341178694174198

vpasolve(EQ,al,30)
ans =
29.051052025226241425009802190818

vpasolve(EQ,al,35)
ans =
35.335792082102904999018466450802

vpasolve(EQ,al,40)
ans =
38.477970385522512262745808241056

In fact, simply looking for zero crossings, we should find roots near each of the following locations:

x = 1:100;
y = double(subs(EQ,al,x));
xl = find(y(1:99).*y(2:100) < 0)';

There will be infinitely many roots, but they will be bracketed as follows:

[xl,xl+1]
ans =
   3     4
   7     8
  10    11
  13    14
  16    17
  19    20
  22    23
  25    26
  29    30
  32    33
  35    36
  38    39
  41    42
  44    45
  47    48
  51    52
  54    55
  57    58
  60    61
  63    64
  66    67
  69    70
  73    74
  76    77
  79    80
  82    83
  85    86
  88    89
  91    92
  95    96
  98    99

Of course, there will be infinitely many solutions. In fact, we can see a pattern:

fzero(@(al) eq(al,k),[3 4])/pi
ans =
       1.2265

fzero(@(al) eq(al,k),[7 8])/pi
ans =
       2.2383

fzero(@(al) eq(al,k),[10 11])/pi
ans =
        3.242

fzero(@(al) eq(al,k),[13 14])/pi
ans =
       4.2439

So the m'th root (beyond the zero root) can be approximated reasonably well as

al(m) ~ pi*(m+1/4)

for integer m. This would make a very good starting value for any numerical scheme, such as fzero. Be careful, since for large values of m, the cosh and sinh terms will blow up, making any numerical computation in double precision impossible. Symbolic computations using vpasolve will be far more stable.

In fact, we can see that as m grows very large, that the root does approach my suggested approximation.

double(vpasolve(EQ,al,pi*(200+1/4))/pi)
ans =
          200.249873456277
>> double(vpasolve(EQ,al,pi*(250+1/4))/pi)
ans =
          250.249898747801
>> double(vpasolve(EQ,al,pi*(500+1/4))/pi)
ans =
          500.249949356665
>> double(vpasolve(EQ,al,pi*(1000+1/4))/pi)
ans =
          1000.24997467402
>> double(vpasolve(EQ,al,pi*(10000+1/4))/pi)
ans =
           10000.249997467

A little mathematical effort would probably show this conclusion holds as m goes to infinity. I have to leave something to the student to prove though.

[Bastian Harvey]

Hi John,

Thank you very much this has been a massive help, just another question, as I am only wanting to find the first 5 eigenvalues so shall i just use:

x = 1:100;

y = double(subs(EQ,al,x));

xl = find(y(1:99).*y(2:100) < 0)';

to find these locations of the roots and then use the vpasolve to find what the eigenvalue is at each location? Sorry if this is very basic stuff, I have just been told to use Matlab for my dissertation so there is a steep learning curve!

Many Thanks,

Bastian

[John D'Errico]

Please don't post an answer just to add a comment or a further question.

Anyway, what you wrote will locate brackets for the first 30 roots or so. As I showed, you can approximate the nth root very well as

pi*(m + 1/4 - p(m))

where p(m) goes to zero. Some mathematical analysis should allow you to show that p(m) indeed goes to zero as m goes to infinity. The exact dependence is not certain to me at this point, but it is your thesis, not mine. :)

So I would simply use

pi*(m + 1/4)

as a starting value. If you want a bracket to insure convergence of a tool like fzero, this would probably suffice:

pi*(m + 1/4) + [-.5, 0.5]

thus,

k = 0.5;
nroots = 5;
sol = zeros(1,nroots);
for m = 1:nroots
  sol(m) = fzero(@(al) eq(al,k),pi*(m + 0.25) + [-0.5 0.5]);
end

sol
ans =
          3.85305048361731          7.03182813669153          10.1850152226347          13.3326619299354          16.4779571341451

sol/pi
ans =
        1.22646406090063          2.23830041385426          3.24199103629701          4.24391810144468           5.2450966599111

[Bastian Harvey]

Sorry I put my comment in the wrong section before but i think i deleted it and put it into the correct place. What is the purpose of dividing through by pi? Is it just to be able to easily approximate a chosen root? Are the eigenvalues not just the roots of the graph without dividing through by pi?

Thank you again! :)

[John D'Errico]

Thank you form moving it. I often move answers into comments, so you saved me that effort.

Dividing by pi allows you to recognize that the solutions are all roughly periodic, of the form that I showed, so approximately non-integer multiples of pi. The eigenvalues themselves are of course not divided by pi, but they are of the form

pi*(m + 1/4 - p(m))

where p(m) is small, and goes to zero for large m. You can choose to use this form for starting values for the roots in a numerical solver, or you can see if mathematical analysis can provide more information about exactly how p(m) goes to zero, so yielding a better approximation yet to the roots.

For example, suppose you decide to expand this function as a Taylor series around the point

    pi*(m + 1/4)

for any general m? Thus it will be an expansion in the form of f(x-pi*(m + 1/4)). Now analyze the behavior of that relation. This would be my first idea in terms of how to show that p(m) goes to zero, as well as predicting the actual behavior to yield a higher order approximation for the m'th root. Of course, this may be more mathematical depth than you want to put into an engineering thesis. (Graduate level engineering work is often just a disguise for applied mathematics though. That said by an applied mathematician who also has a MechE degree.) In the end, it would be a decision you and your thesis advisor may choose to make together.