Solve Non Linear Equations Numerically - Python

Question

I am looking to reproduce results from a research article.

I am at a point, where I have to find the maximum value of the following equation (w), and corresponding independent variable value (k). k is my only variable.

from sympy import *
import numpy as np
import math

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

k = Symbol('k', real=True)
t = 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 + sqrt(w2 + w3 + w4*w5))/w6;

I was able to solve this using Sympy - diff & solve functions, only when I give t = 1 or a any constant.

Do anyone have suggestions on finding the maximum value of this function? Numerically also works - however, I am not sure about the initial guess value. Good thing is I only have one independent variable.

Edit:

As per the answers given here regarding gradient descent, and also plotting and seeing the maximum value. I literally copied the code lines, that include plotting and I got a different plot.

Any thoughts on why this is happening? I am using Python 3.7

Answer 1

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])

Answer 2

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?

Answer 3

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()

Answer 4

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