# How iteratively/numerically solve for a value

**URL:** <https://discourse.julialang.org/t/how-iteratively-numerically-solve-for-a-value/75019>\
**Category:** New to Julia\
**Tags:** roots\
**Created:** [January 22, 2022, 1:43am UTC](https://discourse.julialang.org/t/how-iteratively-numerically-solve-for-a-value/75019 "2022-01-22T01:43:08Z")\
**Posts on this page:** 10\
**Page:** 1

<div class="post-metadata">

**Author:** ![CentraCep](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/centracep/32/32454_2.png) [@CentraCep](https://discourse.julialang.org/u/CentraCep)\
**Post date:** [January 22, 2022, 1:43am UTC](https://discourse.julialang.org/t/how-iteratively-numerically-solve-for-a-value/75019/1 "2022-01-22T01:43:08Z")

</div>

I wrote a program intended to calculate the plasma beta of a hydrostatic solar atmosphere:

```julia
B_0 = 0.02 # base magnetic field [T]
p_0 = 0.015 # base pressure [J/m^-3]
h_D = 7.5E7 # dipole depth [m]
R_☉ = 6.96E8 # solar radius [m]
M_☉ = 1.9891E30 # solar mass [kg]
G = 6.673E-11 # gravitational constant [m^3 kg^-1 s^-2]
μ = 0.61 # mean molecular weight of solar corona
m_H = 1.6726E-27 # mass of H particle [kg]
μ_0 = (4*pi)*10^(-7) # permeability [H/m]
k = 1.3806E-23 # boltzmann constant [J/K]
T = 1E6 # temperature [K]

g = (G*M_☉)/(R_☉^2)

lambda_p = (k*T)/(μ*m_H*g) # coronal scale height for 1MK [m]

p(h) = p_0*exp(-h/lambda_p)
B(h) = B_0*(1 + (h/h_D))^(-3)

```

I am not very familiar with using iterative/numerical methods (let alone implement them with Julia) and was wondering if there was a way to determine the value of `h` when `beta = 1`?

---

<div class="post-metadata">

**Author:** ![cvanaret](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cvanaret/32/11594_2.png) [@cvanaret](https://discourse.julialang.org/u/cvanaret)\
**Post date:** [January 22, 2022, 1:50am UTC](https://discourse.julialang.org/t/how-iteratively-numerically-solve-for-a-value/75019/2 "2022-01-22T01:50:41Z")

</div>

You can use the classical [Newton method](https://en.wikipedia.org/wiki/Newton%27s_method).

---

<div class="post-metadata">

**Author:** ![CentraCep](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/centracep/32/32454_2.png) [@CentraCep](https://discourse.julialang.org/u/CentraCep)\
**Post date:** [January 22, 2022, 3:55am UTC](https://discourse.julialang.org/t/how-iteratively-numerically-solve-for-a-value/75019/3 "2022-01-22T03:55:04Z")

</div>

I was wondering specifically if there was a Julia-esque way, or a way that uses a certain Julia package? Otherwise, I might as well learn Newton method.

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [January 22, 2022, 3:58am UTC](https://discourse.julialang.org/t/how-iteratively-numerically-solve-for-a-value/75019/4 "2022-01-22T03:58:03Z")

</div>

The simplest answer is probably to formulate this as a root finding problem (ie find the root of `f(h) = p(h)/(B(h)^2 / (2*μ_0) )-beta` and use `Roots.jl`

---

<div class="post-metadata">

**Author:** ![CentraCep](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/centracep/32/32454_2.png) [@CentraCep](https://discourse.julialang.org/u/CentraCep)\
**Post date:** [January 22, 2022, 4:32am UTC](https://discourse.julialang.org/t/how-iteratively-numerically-solve-for-a-value/75019/5 "2022-01-22T04:32:38Z")

</div>

So, something like this? I am not familiar with `find_zero`, specifically the second argument.

```julia
f = h -> p(h)/(B(h)^2 / (2*μ_0) ) - beta(h)
Roots.find_zero(f,2*h_D) 

```

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [January 22, 2022, 4:35am UTC](https://discourse.julialang.org/t/how-iteratively-numerically-solve-for-a-value/75019/6 "2022-01-22T04:35:07Z")

</div>

The second argument is an initial guess. You can pass anything you want, but if you happen to have an idea of what the root should be, it can significantly improve performance.

---

<div class="post-metadata">

**Author:** ![CentraCep](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/centracep/32/32454_2.png) [@CentraCep](https://discourse.julialang.org/u/CentraCep)\
**Post date:** [January 22, 2022, 4:40am UTC](https://discourse.julialang.org/t/how-iteratively-numerically-solve-for-a-value/75019/7 "2022-01-22T04:40:52Z")

</div>

Understood, thanks! Anyway, is this an appropriate approach for determining h(beta = 1)? And do I make the result print more decimals?

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [January 22, 2022, 4:45am UTC](https://discourse.julialang.org/t/how-iteratively-numerically-solve-for-a-value/75019/8 "2022-01-22T04:45:46Z")

</div>

it should print all the decimals it can compute. You should be able to get more by doing your calculation with `BigFloat`s (but it will be ~50x slower).

---

<div class="post-metadata">

**Author:** ![Seif\_Shebl](https://avatars.discourse-cdn.com/v4/letter/s/eada6e/32.png) [@Seif\_Shebl](https://discourse.julialang.org/u/Seif_Shebl)\
**Post date:** [January 23, 2022, 11:36pm UTC](https://discourse.julialang.org/t/how-iteratively-numerically-solve-for-a-value/75019/9 "2022-01-23T23:36:20Z")

</div>

Being an ex-MATLAB et al., I really miss the simplicity of doing such things symbolically. Even if a symbolic solution was not found, the algorithm will automatically fallback to a numerical solution without exposing the user to advanced details of the solution algorithm used.

```julia
>> h = solve(2*mu0*p/B^2 - beta)
Warning: Unable to solve symbolically. Returning a numeric solution using vpasolve. 
h =
-234392042.24591558379314753229238

```

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [January 23, 2022, 11:40pm UTC](https://discourse.julialang.org/t/how-iteratively-numerically-solve-for-a-value/75019/10 "2022-01-23T23:40:02Z")

</div>

If you want a symbolic solution, you can try `Symbolics.jl`. I’m not sure how it works though.