# Help for Solving Two Equations For Two Unknows (from a Paper)

I am trying to reproduce numerically the results found in this paper: https://ieeexplore.ieee.org/document/1707778

You don't necessarily need to read it. Basically I boils down to solving a two equations for two unknowns, deflection and voltage: ζ and V

$$\begin{equation*} (k_{0}+k_{1})\zeta+k_{2}\zeta^{3}=\mu f_{e}(\zeta) \tag{16} \end{equation*}$$

\begin{align*} f_{e}(\zeta)((k_{0}+k_{1})&+3k_{2}\zeta^{2})\\ &- \frac{df_{e}(\zeta)}{d\zeta}((k_{0}+k_{1})\zeta+k_{2}\zeta^{3})=0. \tag{18} \end{align*}

By solving (18) for ζ, and by substituting the value of ζ into (16), we determine the deflection and voltage.

Where: $$\begin{equation*} \mu=\frac{6\ell^{4}\epsilon_{0}\epsilon_{r}V^{2}}{Eh^{3}g_{0}^{3}}. \end{equation*}$$

$$\begin{equation*} k_{0}=\int_{0}^{1}(\bar w^{\prime \prime}(x))^{2}dx\\ k_{1}=N_{0} \int_{0}^{1}(\bar{w}^{\prime}(x))^{2}dx\\ k_{2}=N_{1}\left(\int_{0}^{1}(\bar w^{\prime}(x))^{2}dx\right)^{2}\\ f_{e}(\zeta)=\frac{1}{\mu}\int_{0}^{1}F_{e}(\zeta \bar w(x))\bar w(x)dx.\\ \end{equation*}$$

and

$$\begin{equation*} N_{1}=6\frac{g_{0}^{2}}{h^{2}}\\ N_{0}=\frac{N_{0}\ell^{2}}{EI}=12\frac{\tilde{\sigma}\ell^{2}}{Eh^{2}} \end{equation*}$$

$$\begin{equation*} F_{e}=\frac{\mu}{(1-\hat w)^{2}}\left(\mathcal{F}+\eta\frac{\partial \mathcal{F}}{\partial\eta}\right) \end{equation*}$$

$$\begin{equation*} \mathcal{F}(\beta,\ \eta)=1-0.36\frac{\beta}{\eta}+0.85(\frac{\beta}{\eta})^{0.76}+2.5\frac{\beta}{\eta^{0.76}}. \end{equation*}$$

where: $$\begin{equation*} \beta=\frac{h}{b},\ \eta=\frac{h}{g}. \tag{6} \end{equation*}$$

I am trying to implement it into Mathematica, but it is not working.

eq16[ζ_]:= -μf[ζ]+ζ^3 k2+ζ (k0+k1)
eq18[ζ_]:= f[ζ]((k0+k1)+3 k2 ζ^2)-D[f[ζ],ζ]((k0+k1)ζ+k2 ζ^3)

ℱ[η_,β_]:= 1 - 0.36 β/η + 0.85 (β/η)^0.76 + 2.5 β/η^0.76
F[w_]:= μ/(1-w)^2(ℱ[η,β]+η*D[ℱ[η,β],η])

f[ζ_] := 1/μ Integrate[F[ζ*w[x]]*w[x],{x,0,1}]

μ = (6*L^4*eps0*epsR*V^2)/(e*h^3*g^3)
k0 = Integrate[D[w[x],{x,2}]^2,{x,0,1}]
k1 = n0*Integrate[D[w[x],{x,1}]^2,{x,0,1}]
k2 = n1*(Integrate[D[w[x],{x,1}]^2,{x,0,1}])^2

n0 = 12 (sig*L^2)/(e h^2)
n1 = 6 g^2/h^2

(*Numeric Parameter*)
L = 100*10^-6;
b = 1*10^-6;
h = 2*10^-6;
g = 4*10^-6;
e = 169*10^9;
eps0 = 8.85*10^-12
epsR = 1
sig = 100*10^6*(1-0.066)
β=h/b
η=h/g

(*Trial Function*)
w[x] := 16 x^2 (1 - x)^2;

(*Solving the System*)
ζSolve = ζ/.Solve[eq18[ζ]==0,ζ]
VSolve = V/.Solve[eq16[ζSolve]==0,V]


It just keeps running and running and running... I guess that I need to solve it somehow numerically. In fact in the paper they write:

We emphasize that, once the trial function w(x) has been chosen, (18) reduces to a nonlinear algebraic equation for ζ, where the derivative df(ζ)/dζ can be computed for any ζ by numerical integration. In general, (16) cannot be solved analytically, and standard root finding techniques such as the bisection algorithm may be applied

...but how do you do that in Mathematica ?

Any help would be very much appreciated !!

P.S. The result for V should be around: 945

• You have constants that differ by dozens of orders of magnitude. Numerical results will likely be very bad. Please, work in natural units, where all parameters are of order 1. – AccidentalFourierTransform Dec 3 '18 at 19:58
• @AccidentalFourierTransform This is true, but the author of the paper already thought about it, and hence all equations are dimensionless, so it should not be a problem. – james Dec 3 '18 at 20:08
• In your definition of eq16, do you want a minus sign in front of the $\mu f$ term? – LouisB Dec 3 '18 at 21:00
• @LouisB Very true. Thanks for the catch ! – james Dec 3 '18 at 21:21
• @james parameters not defined η,β. what are they needed for? – Alex Trounev Dec 3 '18 at 23:13

## 2 Answers

I read the article, made all the corrections.I managed to find the correct model that reproduces the data from table 3 of Romesh C. Batra, ASME, Fellow, Maurizio Porﬁri, and Davide Spinello. Electromechanical Model of Electrically Actuated Narrow Microbeams. JOURNAL OF MICROELECTROMECHANICAL SYSTEMS, VOL. 15, NO. 5, OCTOBER 2006, p. 1175-1189

    eq16[\[Zeta]_,
V_] := -\[Mu][V]*f[\[Zeta]] + \[Zeta]^3 k2 + \[Zeta] (k0 + k1)
eq18[\[Zeta]_] :=
f[\[Zeta]] ((k0 + k1) + 3 k2 \[Zeta]^2) -
D[f[\[Zeta]], \[Zeta]] ((k0 + k1) \[Zeta] + k2 \[Zeta]^3)
\[Mu][V_] := (6*L^4*eps0*epsR*V^2)/(e*h^3*g0^3)
\[ScriptCapitalF][\[Eta]_, \[Beta]_] :=
1 - 0.36 \[Beta]/\[Eta] + 0.85 (\[Beta]/\[Eta])^0.76 +
2.5 \[Beta]/\[Eta]^0.76
F[x_] := 1/(1 - x)^2

(*Trial Function*)
w[x_] := 16 x^2 (1 - x)^2;

(*Numeric Parameter*)
L = 100*10^-6;
b = 1*10^-6;
h = 2*10^-6;
g0 = 4*10^-6;
e = 169*10^9;
eps0 = 8.854187817*10^-12;
epsR = 1;
nu = .066;
sig0 = 100*10^6;
sig = (1 - nu)*sig0;
p = {n0 = 12 (sig*L^2)/(e h^2),
n1 = 6 g0^2/h^2, k0 = NIntegrate[D[w[x], {x, 2}]^2, {x, 0, 1}],
k1 = n0*NIntegrate[D[w[x], {x, 1}]^2, {x, 0, 1}],
k2 = n1*(NIntegrate[D[w[x], {x, 1}]^2, {x, 0, 1}])^2};
B = (\[ScriptCapitalF][\[Eta], \[Beta]] + \[Eta]*
D[\[ScriptCapitalF][\[Eta], \[Beta]], \[Eta]]) /. {\[Eta] ->
h/(g0*(1 - \[Zeta]*w[x])), \[Beta] -> h/b} // Simplify;
G = Interpolation[
Table[{\[Zeta],
NIntegrate[
B*F[\[Zeta]*w[x]]*w[x], {x, 0, 1}]}, {\[Zeta], -0.5, .75, .01}]]
f[\[Zeta]_] := G[\[Zeta]]

(*Solving the System*)
\[Zeta]0 = FindRoot[eq18[\[Zeta]] == 0, {\[Zeta], .62}]

(*Out[]= {\[Zeta] -> 0.620671}*)

FindRoot[eq16[\[Zeta], V] /. \[Zeta]0, {V, 1}]

(*Out[]= {V -> 945.44}*)

• Thanks a lot for your answer. Great idea with the interpolation ! Unfortunatly, when I run it, it does the interpolation, but somehow FindRoot does not work. My output is: "FindRoot[eq18[[Zeta]] == 0, {[Zeta], 0.4}]". It does not seem to evaluate it. – james Dec 4 '18 at 7:58
• I see that you put some constants (n0,n1,k1,k2) into brackets. Why ? – james Dec 4 '18 at 8:13
• @james I put the parameters in brackets for easier typing while debugging. What is your version, operating system and machine? – Alex Trounev Dec 4 '18 at 11:38
• I have Windows 10 and Mathematica 11.3 – james Dec 5 '18 at 5:57
• I also checked the article and I come to the same conclusion as you regarding the typo. I will try to check the math in the paper to find the error... but this could take a while – james Dec 5 '18 at 5:59

Udated numerical data:

Clear["Global*"]

(*Trial Function*)
w[x] := 16 x^2 (1 - x)^2

eq16[ζ_] := -μ f[ζ] + ζ^3 k2 + ζ (k0 + k1)
eq18[ζ_] := f[ζ] ((k0 + k1) + 3 k2 ζ^2) - D[f[ζ], ζ] ((k0 + k1) ζ + k2 ζ^3)

ℱ[η_, β_] := 1 - 0.36 β/η + 0.85 (β/η)^0.76 + 2.5 β/η^0.76
F[ww_] := μ/(1 - ww)^2 (ℱ[η, β] + η*D[ℱ[η, β], η])


I used ww since we already have a w[x]. Now what took forever for me was the following integration for x from 0 to 1. Indefinite integration works much faster and then manually apply the limits.

fz[x_] = 1/μ Integrate[F[ζ*w[x]]*w[x], x] // Simplify
(*Long answer in RootSum objects*)

f[ζ_] := fz - fz // Simplify

μ = (6*L^4*eps0*epsR*V^2)/(e*h^3*g^3)
k0 = Integrate[D[w[x], {x, 2}]^2, {x, 0, 1}]
k1 = n0*Integrate[D[w[x], {x, 1}]^2, {x, 0, 1}]
k2 = n1*(Integrate[D[w[x], {x, 1}]^2, {x, 0, 1}])^2

n0 = 12 (sig*L^2)/(e h^2)
n1 = 6 g^2/h^2

(*Numeric Parameter*)
L = 100*10^-6;
b = 1*10^-6;
h = 2*10^-6;
g = 4*10^-6;
e = 169*10^9;
eps0 = 8.85*10^-12
epsR = 1
sig = 100*10^6*(1 - 0.066)
β = h/b
η = h/g


Find ζ

Plot[Evaluate[eq18[ζ] // N // Chop], {ζ, -.5, .7}, PlotRange -> All] FindRoot[Re[eq18[ζ]] == 0, {ζ, .55}] // Chop
(*{ζ -> 0.531648}*)

ζ = ζ /. %

Solve[(eq16[ζ] // Chop) == 0, V]
(*{{V -> -795.934}, {V -> 795.934}}*)
`
• Very nice answer ! Thanks a lot ! I updated the numeric constants since I did a mistake when copying from the paper. However, the result is still not the same... – james Dec 4 '18 at 8:14
• @james What should be the result? – Alex Trounev Dec 4 '18 at 11:46
• @AlexTrounev It should be 945. I also spotted a little copy mistake in ℱ[η_,β_] but it is now fixed in the question. – james Dec 4 '18 at 13:22
• @james then we have to check all the data more carefully. – Alex Trounev Dec 4 '18 at 14:07
• Hi ! So I have contacted the author of the paper. He told me that we have to adapt η because the gap decreases with increasing voltage. Basically on needs to replace η = h/g with η = h/(g-w[x]). I tried it with your code, but I don't get any result. – james Dec 6 '18 at 10:59