# eigenvalue and eigenfunction of complex euler-bernoulli beam

So I am new to Mathematica and am trying to solve the euler-bernoulli modal equation for a U-shaped Cantilever beam given by equations :-

where i is the index of the region. In total there are 2 regions, each with its own EI and mu values respectively. Region 1 spans from x = 0 to x = Lleg and region 2 spans from x = Lleg to x = L. The solution is given by the expression :-

and the boundary conditions are as follows :-

I know mathematica has NDEigensystem function which can help me with this but I don't know how to use it correctly.

Edit :- I would also Like to develop an analytical expression of Phi(x) as a function of x for the 2 regions since I need to integrate that expression to obtain some discrete parameters as follows :-

The code block is as follows :-

EAu = 78*10^9; (*Youngs Modulus of Gold*)
ESiN = 250*10^9; (*Youngs Modulus of Silicon Nitride*)
rhoAu = 19300; (*Density of Gold*)
rhoSiN = 3440; (*Density of Silicon Nitride*)
b11 =1.5; (*width of gold, section I*)
b12 = 4.5; (*width of gold, section II*)
b21 = b11; (*width of SiN, section I*)
b22 = b12; (*width of SiN, section II*)
h11 = 20*10^(-3); (*height of gold, section I*)
h21 = 510*10^(-3); (*height of SiN, section I*)
h12 = h11; (*height of gold, section II*)
h22 = h21; (*height of SiN, section II*)
IAu1 =(1/12)*b11*h11^3; (*2nd Moment of Area, gold, section I, about   the center*)
IAu2 = (1/12)*b12*h12^3; (*2nd Moment of Area, gold, section II, about the center*)
ISiN1= (1/12)*b21*h21^3; (*2nd Moment of Area, SiN, section I, about the center*)
ISiN2 = (1/12)*b22*h22^3; (*2nd Moment of Area, SiN, section II, about the center*)

EIsys1 = 2*EAu*(IAu1 + b11*h11*(0.5*(h11+h21)-0.5*h11)^2) + 2*ESiN*(ISiN1 + b21*h21*(0.5*(h11+h21)-0.5*h21)^2)
EIsys2 = EAu*(IAu2 + b12*h12*(0.5*(h12+h22)-0.5*h12)^2) + ESiN*(ISiN2 + b22*h22*(0.5*(h12+h22)-0.5*h22)^2)

musys1 = 2*rhoAu*b11*h11 + 2*rhoSiN*b21*h21 (*mass per unit length, section I*)
musys2 = rhoAu*b12*h12 + rhoSiN*b22*h22 (*mass per unit length, section II*)

AR = 5; (*Input Value, Aspect Ratio of Beam*)
L = AR*b12 (*Length of Beam, total*)
Lleg = AR*b11 (*Length of Beam, Section I*)

EIL = EIsys1
EIR = EIsys2
\[Mu]L = musys1
\[Mu]R = musys2
bleg = b11
b = b12
m = Lleg
eqnL = EIL \[Phi]L''''[x] - \[Mu]L *(\[Omega]^2)* \[Phi]L[x] == 0
eqnR = EIR \[Phi]R''''[x] - \[Mu]R *(\[Omega]^2)* \[Phi]R[x] == 0

bcs = {\[Phi]L[0] == 0, \[Phi]L'[0] == 0,
\[Phi]L[m] == \[Phi]R[m], \[Phi]L'[m] == \[Phi]R'[m],
2 bleg \[Phi]L''[m] == b \[Phi]R''[m], 2 bleg \[Phi]L'''[m] == b \[Phi]R'''[m],
\[Phi]R''[L] == 0, \[Phi]R'''[L] == 0}

• Maybe you can make a few steps towards solution. First, get rid of times and formulate an eigenvalue equation. Try to type this equation and boundary conditions in Mathematica form to help people who want to help you. Then look at NDEigensystem. – yarchik Oct 8 '20 at 15:31
• Do you want to solve this numerically? – SPPearce Oct 11 '20 at 16:18
• I have a package available for finding eigenvalues numerically, which can handle this kind of interface problem. See mathematica.stackexchange.com/questions/180005/… for an example. If it is interesting to you, I can show how to use it for this problem. – SPPearce Oct 11 '20 at 16:33
• You can build a full finite element model in Mathematica and then find the eigenvalues and vectors. This post gives an example. Why do you want to joint together a set of beams? – Hugh Oct 12 '20 at 11:52
• Please can you post Mathematica code so that we can copy and don't have to type it into a notebook? Just paste the code, select it and click the curly brackets. – Hugh Oct 15 '20 at 16:37

I have a package that implements solving eigenvalue problems, including interface problems like this.

First we need to install (first time only):

Needs["PacletManager"]
PacletInstall["CompoundMatrixMethod",
"Site" -> "http://raw.githubusercontent.com/paclets/Repository/master"]


Needs["CompoundMatrixMethod"]


We convert the system of ODEs into a matrix form via my function ToMatrixSystem:

sys = ToMatrixSystem[{eqnL, eqnR}, bcs, {ϕL, ϕR}, {x, 0, m, L}, ω];


The method generates something called the Evans function, roots of which correspond to eigenvalues of the original system.

This can be evaluated for a given value of $$\omega$$, say $$\omega = 1$$, with:

Evans[1, sys]
(* 4.54519 *)


This is not zero, so $$\omega = 1$$ is not an eigenvalue of this equation. Also note that it doesn't get fooled by $$\omega = 0$$, which the determinant will vanish at.

We therefore just need to find roots of this function, via plotting or FindRoot.

FindRoot[Evans[ω, sys], {ω, 1}]
(* {ω -> 6.79439} *)


And you can see multiple roots in a plot:

Plot[Evans[ω, sys], {ω, 0, 500}]


• Hi, Thanks a lot for the input. I am using this for microstructure so the dimensions are small. I get the following error imgur.com/PYFNCZh which yields incorrect freq. Also, is there a way to extract mode shape function Phi(x) for a given eigfreq? – Saransh Oct 12 '20 at 10:50
• FYI, these are the parameters I used. imgur.com/i5dBTdk – Saransh Oct 12 '20 at 10:58
• @Saransh, you probably want to rescale your units, values of $10^{-15}$ aren't very helpful. Can you edit your question with the actual values you are interested in (also needs a value for $L$). – SPPearce Oct 12 '20 at 12:08
• I have never got round to figuring out how to return the eigenfunction from my package for this method as yet (and I have no time to work on this at the moment). – SPPearce Oct 12 '20 at 12:20
• That is to say I don't have time to work on the general case, but I might be able to get code out for the eigenfunctions for one particular case. – SPPearce Oct 12 '20 at 18:37

parms = {EIL -> 4.31671*10^(-15), EIR -> 1.29501*10^(-14), \[Mu]L -> 3.2106*10^(-9), \[Mu]R -> 9.6318*10^(-9), bleg -> 1.5*10^(-6), b -> 4.5*10^(-6), m -> 7.5*10^(-6), L -> 22.5 10^(-6)};
eqnL = \[Phi]L''''[x] - \[Mu]L /EIL \[Omega]^2 \[Phi]L[x] == 0;
eqnR = \[Phi]R''''[x] - \[Mu]R /EIR  \[Omega]^2 \[Phi]R[x] == 0;
solL = DSolve[eqnL, \[Phi]L, x][[1]];
solR = DSolve[eqnR, \[Phi]R, x][[1]];
\[Phi]Lx = \[Phi]L[x] /. solL;
\[Phi]Rx = \[Phi]R[x] /. solR /. {C[1] -> C[5], C[2] -> C[6], C[3] -> C[7], C[4] -> C[8]};
equ1 = \[Phi]Lx /. {x -> 0};
equ2 = D[\[Phi]Lx, x] /. {x -> 0};
equ3 = (\[Phi]Lx - \[Phi]Rx) /. {x -> m};
equ4 = D[\[Phi]Lx - \[Phi]Rx, x] /. {x -> m};
equ5 = D[2 bleg \[Phi]Lx - b \[Phi]Rx, {x, 2}] /. {x -> m};
equ6 = D[2 bleg \[Phi]Lx - b \[Phi]Rx, {x, 3}] /. {x -> m};
equ7 = D[\[Phi]Rx, {x, 2}] /. {x -> L};
equ8 = D[\[Phi]Rx, {x, 3}] /. {x -> L};
M = Grad[{equ1, equ2, equ3, equ4, equ5, equ6, equ7, equ8}, Table[C[k], {k, 1, 8}]];
det = Det[M] /. parms;


Plotting the graphics for $$\det(\omega)$$ we have

gr0 = LogLogPlot[det, {\[Omega], 0, 10^9}, PlotStyle -> {Thick, Blue}]


from which we obtain the two first characteristic frequencies as follows

r1 = Quiet@FindRoot[det == 0, {\[Omega], 6.3 10^6}];
r1a = Quiet@FindRoot[det == 0, {\[Omega], 10^7 }];
r2 = Quiet@FindRoot[det == 0, {\[Omega], 45 10^6 }];
r2a = Quiet@FindRoot[det == 0, {\[Omega], 5 10^7 }];

omega1 = \[Omega] /. r1
omega1a = \[Omega] /. r1a
omega2 = \[Omega] /. r2
omega2a = \[Omega] /. r2a

• Yes, this works fine for DEs that have an analytic solution. You do need to be careful to make sure that the roots of the determinant that you return are not due to points where the solutions to the DEs become linearly dependent. – SPPearce Oct 16 '20 at 7:33
• Hi @Cesareo, how do I get Phi[x]L and Phi[x]R from this? – Saransh Nov 27 '20 at 8:25