python - 数值求解非线性方程 - Python

Question

我希望重现一篇研究文章的结果。

我现在必须找到以下等式 (w) 的最大值，以及相应的自变量值 (k)。 k 是我唯一的变量。

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;

我能够使用 Sympy - diff & solve 函数解决这个问题，只有当我给出 t = 1 或任何常量时。

有人对寻找这个函数的最大值有什么建议吗？数值上也有效——但是，我不确定初始猜测值。幸好我只有一个自变量。

编辑:

根据此处给出的有关梯度下降的答案，以及绘制和查看最大值。我从字面上复制了代码行，其中包括绘图，但我得到了一个不同的绘图。

有没有想过为什么会这样？我正在使用 Python 3.7

Answer 1

有很多方法可以做到这一点。 scipy 特别有一堆优化算法。我将使用梯度下降(或者更恰本地说，梯度 ascent)和 autograd，因为它可能很有趣。

首先，让我们导入 autograd 并将您的函数转换为可调用函数。

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

现在，我们可以使用autograd 来计算w 相对于k 的梯度。您可以添加一些逻辑，以确保一旦我们达到某个容差阈值，该过程就会终止。

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)

现在，让我们绘制结果以确保找到最大值

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