Problem Description:

In solving the "Shifted Neutral Axis" method for the stress distribution of a vertical steel vessel supported on a concrete foundation, we can derive the following expression:

$$ \frac{\left(\gamma+\alpha\right)-\tan (\alpha)}{\left(\gamma+\alpha\right) \sec (\alpha )-\sin (\alpha )}=-\beta $$


$\alpha$ is the angle of the "shifted" axis of rotation ranging from $0$ to $\pi$

$\beta$ is a ratio of dead to lateral load from $0$ to $1$. Real life values usually range from $0.10$ to $0.50$

$\gamma$ is a positive geometric and material constant. Real life values can range from $0.05$ to $0.20$

Essentially, this expression suggests that as the lateral load grows, the neutral axis shifts from a maximum value of $\alpha=\pi$ towards an asymptotic, smaller non-zero value.


Usually, Engineers are given geometric and load constraints ($\beta$ and $\gamma$) and are told to find the system's response, $\alpha$.

The above expression, being implicit, leads itself to a simple solution via FindRoot[], but it would be ideal to solve for $\alpha$ as an explicit function of $\beta$ and $\gamma$, or at the very least come up with an approximated expression to calculate $\alpha$ directly.

I am not sure where to start solving this problem with Mathematica, but pointers would be appreciated.

Edit: The ultimate objective is to obtain an expression that reasonably approximates the shifted neutral axis angle--one that engineers can easily program into an Excel spreadsheet.


2 Answers 2


Generating a differential equation, you can calculate an interpolating function alphasol depending an beta and gamma.

Since i use version 8.0, i had to make an unusal construction for NDSolve. With higher versions use ParametricNDSolve with beta as parameter.

eq1[α_, β_, γ_] = 
eq = ((γ + α) - 
Tan[α])/((γ + α) Sec[α] - 
Sin[α]) == -β

beta[α_] = -eq[[1]] /. γ -> 1/10

invbeta[β_] = InverseFunction[beta][β]

deqγ = 
   D[eq1[α[β, γ], β, γ], γ] // 

αsol = α /. 
  First@NDSolve[{deqγ, α[β, 1/10] == 
  invbeta[β]}, α, {β, .1, .5}, {γ, .05, \

Get alpha-values directly, e.g. αsol[.1, .1] and compare with ContourPlot.

αsol[.1, .1]

ContourPlot[((γ + α) - 
 Tan[α])/((γ + α) Sec[α] - 
 Sin[α]) == -β, {γ, .05, 0.2}, {α, 0, 
1.5}, AspectRatio -> 1, ImageSize -> 300], {β, .1, .5, 
Appearance -> "Labeled"}], 
Plot[αsol[β, γ], {γ, .05, .2}, 
PlotRange -> {0, 1.5}, Frame -> True, AspectRatio -> 1, 
ImageSize -> 300], {β, .1, .5, Appearance -> "Labeled"}
  • $\begingroup$ This is a neat way of getting plots from the interpolating function. The dream would be to generate an expression (analytically derived or fit) that other engineers can program into a spreadsheet. $\endgroup$
    – felimz
    Commented Jul 6, 2020 at 1:32

This is a transcendental equation and in general they don't have analytic solutions. There are special cases, of course, where a "trick" works. But in general no.

See https://en.wikipedia.org/wiki/Transcendental_equation for example.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.