前言
工作中,遇到工程上的一个四元的非线性方程组需要求解,经过各路大神的协助,终将该问题解决,在此进行记录,同时也写给需要的你们。
一、问题描述
二、求解算法
三、代码实现
import numpy as np
from sympy import symbols, diff
from scipy import linalg
import pandas as pd
def fun(x, points):
"""
求方程组(值表示)在点x处的值
:param x: 数值解点
:param points: 给定点(和值,测试时用)
:return:
"""
f = np.zeros(4, dtype=float)
f[0] = x[0] * (x[1] + points['t1']) ** x[3] - points['u1']
f[1] = x[0] * (x[1] + points['t1'] + x[2]) ** x[3] - points['u2']
f[2] = x[0] * (x[1] + points['t2']) ** x[3] - points['u3']
f[3] = x[0] * (x[1] + points['t2'] + x[2]) ** x[3] - points['u4']
return f
def fun_dfun(x, points):
"""
方程组和jacobi矩阵(符号表示)及其在点x处的值
:param x: 分量分别是未知数的一组值
:param points: 给定点(和值,测试时用)
:return: 方程组
"""
# 未知数符号
G = symbols('G')
R = symbols('R')
N = symbols('N')
alpha = symbols('alpha')
# 非线性方程组(符号表示)
f = []
f.append(G * (R + points['t1']) ** alpha - points['u1'])
f.append(G * (R + points['t1'] + N) ** alpha - points['u2'])
f.append(G * (R + points['t2']) ** alpha - points['u3'])
f.append(G * (R + points['t2'] + N) ** alpha - points['u4'])
# 求Jacobi矩阵(符号表示形式)
df = []
for i in range(4):
df.append([diff(f[i], G), diff(f[i], R), diff(f[i], N), diff(f[i], alpha)])
f_x, df_x = np.zeros(4, dtype=float), np.zeros((4, 4), dtype=float)
# 方程组和Jacobi矩阵的值
for i in range(4):
f_x[i] = float(f[i].evalf(subs={G: x[0], R: x[1], N: x[2], alpha: x[3]})) # 值代入方程中
for i in range(4):
for j in range(4):
df_x[i][j] = float(df[i][j].evalf(subs={G: x[0], R: x[1], N: x[2], alpha: x[3]})) # 值带入Jacobi矩阵中
# df_1x = np.linalg.inv(np.array(df_x)) # Jacobi矩阵得逆, 太耗时
return f_x, df_x
def newton(x, points):
"""
用牛顿迭代法求解最优值
:param x: 初始值
:param points: 给定点(和值,测试时用)
:return: 最优值
"""
count = 1
while count < 1000:
f_x, df_x = fun_dfun(x, points)
delta_x = linalg.solve(df_x, -f_x) # 求解线性方程组
x1 = x + delta_x
error = np.sum(np.square(fun(x1, points))) # 迭代中两次值得欧氏距离
if error > 1e-8:
x = x1
count += 1
else:
break
if count < 1000:
print("迭代次数为: ", count)
print("G的最优解为: ", x[0])
print("R的最优解为: ", x[1])
print("N的最优解为: ", x[2])
print("alpha的最优解为:", x[3])
else:
print("迭代失败")
return x
四、测试
def main():
# 测试用例
t = [1, 2, 4, 5]
u = [2 * (4 + i) ** 0.7 for i in t] # G: 2, N: 1, R:3, alpha:0.7
points = {}
points['t1'] = t[0]
points['t2'] = t[2]
points['u1'] = u[0]
points['u2'] = u[1]
points['u3'] = u[2]
points['u4'] = u[3]
x = np.ones(4, dtype=float) # 初值
a = newton(x, points)
print(a)
测试结果如下:
END
参考资料
-(需要补)