Solve Non Linear Equations Numerically - Python

There are a bunch of ways to do this. scipy in particular has a bunch of optimization algorithms. I'm going to use gradient descent (or, perhaps more appropriately, gradient ascent) and autograd because it might be fun.

First, let's import autograd and turn your function into a callable function.

import autograd.numpy as anp
from autograd import grad
import matplotlib.pyplot as plt


def w(k):
    rho_l = 1352;
    rho_g= 1.225;
    sigma = 0.029;
    nu = 0.02;
    Q = rho_g/ rho_l;
    u = 99.67;
    h = 1.6e-5; # Half sheet thickness
    t = anp.tanh(k*h);
    w1 = -2 * nu * k**2 * t ;
    w2 = 4* nu**2 * k**4 * t**2;
    w3 = - Q**2 * u**2 * k**2;
    w4 = -(t + Q)
    w5 = (- Q* u**2 * k**2  +  (sigma*k**3/ rho_l));
    w6 = -w4;

    w = ( w1 + anp.sqrt(w2 + w3 + w4*w5))/w6;

    return w

Now, we can use autograd to compute the gradient of w with respect to k. You can add some logic to ensure that the procedure will terminate once we meet some tolerance threshold.

dwdk = grad(w)

#Do gradient descent

k = 10.0 #An initial guess
learning_rate = 1
for i in range(1000):

    k+=learning_rate*dwdk(k)

And now, let's plot the result to ensure we found the maximum

K = np.arange(0,1000)
plt.plot(K,w(K))
plt.scatter(k, w(k), color = 'red')
plt.show()

My thought is that the function returning values for w (your blue line) might be a truncated estimate (with a polynomial maybe)? Is the formula for w off by a factor of 10 or so?

Here is another method. It is an implementation of the Metropolis algorithm, a so-called Markov Chain Monte Carlo method. Using the definition of w, it is possible to construct a Markov Chain of w(k) called wlist. The tail of this chain should be the maximum of w, and we can recover k that got it by storing the k values in a list called kvalues.

import math
import random

klist = [1.0]
wlist = [w(1.0)] # initialize the chain

# you can tune the value passed to `range`
for _ in range(5000):

    k = random.gauss(klist[-1], 0.2*klist[-1]) # q

    if k <= 0.0: # assuming max has positive `k` arg
        continue

    w_hat = w(k)

    if w_hat > wlist[-1]:
        klist.append(k)
        wlist.append(w_hat)

    else:
        u = random.random()
        try:
            alpha = math.exp(-w_hat) / math.exp(-wlist[-1])
        except ZeroDivisionError:
            alpha = 1.0
        if u >= alpha:
            klist.append(k)
            wlist.append(w_hat)
        else:
            klist.append(klist[-1])
            wlist.append(wlist[-1])

wlist[-10:], klist[-10:]

Which should return approximately (my seed is not set) something like this:

([8679.594992731532,
  8679.594992731532,
  8679.594992731532,
  8679.594992731532,
  8679.594992731532,
  8679.594992731532,
  8679.594992731532,
  8679.594992731532,
  8679.594992731532,
  8679.594992731532],
 [416.22335719432436,
  416.22335719432436,
  416.22335719432436,
  416.22335719432436,
  416.22335719432436,
  416.22335719432436,
  416.22335719432436,
  416.22335719432436,
  416.22335719432436,
  416.22335719432436])

I do not expect that there is an analytical solution to this problem. There is a theory of pfaffian functions for which it is possible to provide certificate that there are no roots for given range. See https://en.m.wikipedia.org/wiki/Pfaffian_function. However this heavy artilery.

If you are not sure about initial guess try to compute the function for a large number of random points (e.g.million) and select the best one as starting point. This approach works really well for low dimensional, differentiable problems