# 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. Dec 3, 2018 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. Dec 3, 2018 at 20:08
• In your definition of eq16, do you want a minus sign in front of the $\mu f$ term? Dec 3, 2018 at 21:00
• @LouisB Very true. Thanks for the catch ! Dec 3, 2018 at 21:21
• @james parameters not defined η,β. what are they needed for? Dec 3, 2018 at 23:13

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. Dec 4, 2018 at 7:58
• I see that you put some constants (n0,n1,k1,k2) into brackets. Why ? Dec 4, 2018 at 8:13
• @james I put the parameters in brackets for easier typing while debugging. What is your version, operating system and machine? Dec 4, 2018 at 11:38
• I have Windows 10 and Mathematica 11.3 Dec 5, 2018 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 Dec 5, 2018 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

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... Dec 4, 2018 at 8:14
• @james What should be the result? Dec 4, 2018 at 11:46
• @AlexTrounev It should be 945. I also spotted a little copy mistake in ℱ[η_,β_] but it is now fixed in the question. Dec 4, 2018 at 13:22
• @james then we have to check all the data more carefully. Dec 4, 2018 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. Dec 6, 2018 at 10:59