# 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 '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. Dec 3 '18 at 20:08
• In your definition of eq16, do you want a minus sign in front of the $\mu f$ term? Dec 3 '18 at 21:00
• @LouisB Very true. Thanks for the catch ! Dec 3 '18 at 21:21
• @james parameters not defined η,β. what are they needed for? Dec 3 '18 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 '18 at 7:58
• I see that you put some constants (n0,n1,k1,k2) into brackets. Why ? 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? Dec 4 '18 at 11:38
• I have Windows 10 and Mathematica 11.3 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 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

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 '18 at 8:14
• @james What should be the result? 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. Dec 4 '18 at 13:22
• @james then we have to check all the data more carefully. 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. Dec 6 '18 at 10:59