如何用Python 和牛顿法解四元一次方程组

Posted 2023-03-19

我是试图用Sympy这个库来解决问题，却出错了。
from math import *
from pylab import *

#Avoid confict betwwen sympy and math
import sympy as spy

from WoyTools import *

r1=300
r2=1200
r3=800
r31=476.74
r32=1000
r4=800
r5=400
r6=1746.42
r7=600
theta1=pi/2
theta7=pi/2

alfa=float(float(100.02)/180.0*spy.pi)
gama=float(float(28)/180.0*spy.pi)
theta6=float(float(66.371)/180.0*spy.pi)

#set angle velocity of theta2
w1 = float(2*spy.pi)

#set angle acceleration of theta2
a1 = float(0/180.0*spy.pi)

ri=float(180/180.0*spy.pi)

theta5=float(323.13/180.0*spy.pi)

Ax=0
Ay=0
Ex=70
Ey=160
Fx=130
Fy=160

#guess value
theta2=float(float(84.982)/180.0*spy.pi)
theta3=float(float(22.676)/180.0*spy.pi)
theta4=float(float(159.626)/180.0*spy.pi)
theta5=float(float(152.334)/180.0*spy.pi)

time = 1
count=0
for t in arange(0,time+0.1,0.01):
theta1 = w1 * t + 0.5 * a1 * t * t

theta22=spy.Symbol('theta22')
theta33=spy.Symbol('theta33')
theta44=spy.Symbol('theta44')
theta55=spy.Symbol('theta55')

f1=r1*cos(theta1)+r2*spy.cos(theta22)+r31*spy.cos(theta33)*spy.cos(alfa)-r31*spy.sin(theta33)*spy.sin(alfa)+r4*cos(theta44)-r6*cos(theta6)
f2=r1*sin(theta1)+r2*spy.sin(theta22)+r31*spy.sin(theta33)*spy.cos(alfa)+r31*spy.cos(theta33)*spy.sin(alfa)+r4*sin(theta44)-r6*sin(theta6)

f3=r32*spy.cos(theta33)*spy.cos(gama)+r32*spy.sin(theta33)*spy.sin(gama)-r4*cos(theta44)+r5*cos(theta55)+r7
f4=r32*spy.sin(theta33)*spy.cos(gama)-r32*spy.cos(theta33)*spy.sin(gama)-r4*sin(theta44)+r5*sin(theta55)-r6*sin(theta6)

theta2,theta3,theta4,theta5 = spy.nsolve((f1,f2,f3,f4),(theta22,theta33,theta44,theta55),(theta2,theta3,theta4,theta5))

错误提示是在f1那里raise ValueError("Symbolic value, can't compute")
ValueError: Symbolic value, can't compute

有没有其他能用的解牛顿法的库或者如果用Sympy解决这个问题？
Sympy这个库里关于使用牛顿法的语法规则就是要加引号的，我这么写在其他程序里没有问题

比较弱的问一下，你确定不是
'''
theta22=spy.Symbol('theta22')
theta33=spy.Symbol('theta33')
theta44=spy.Symbol('theta44')
theta55=spy.Symbol('theta55')
'''
这段有问题？
多了引号？或者。。 参考技术A 错误的意思不就是：不能计算符号值吗？同意楼上的判断

Python | 牛顿法解一元方程

本文首发于微信公众号："算法与编程之美"，欢迎关注，及时了解更多此系列文章。

问题描述

在日常生活中面对一元一次方程、一元二次方程甚至是一元三次方程时，我们都可以应用所学的数学知识比如因式分解和求根公式去将其轻易地解开。但并不是所有的方程都能进行因式分解或有求根公式又或者是求根公式复杂，这些问题都导致我们解方程时求解困难，形如x-cos(x)=0、e^x * sin(x)-cos(x)=0、x^4–2x^2+x=0等超越方程和高次方程我们人工计算就无法求得其解。

解决方案

我们可以应用牛顿法解决上述问题。牛顿法又称切线法，简单来说就是不断求函数图像的切线与x轴的交点，来逐渐接近函数解的一个迭代过程，其核心思想就是不断地逼近函数的解。

牛顿法的步骤如下：

.选取一个初始迭代点，如‘1’；

.根据公式求出下一个接近零点的点；

③重复步骤②直到满足：

1). (eps为人设的精度要求)

2).达到最大迭代次数

牛顿法的原理如图：

通过上图我们可以发现，每迭代一次的值都将更接近方程的解。

接下来我们用Python代码来实现牛顿法并求方程的解：

from sympy import * # 定义方程 def function():     x = symbols('x') # 符号变量的定义     fx = -exp(x) * sin(x) + cos(x) - 1 return fx   # 求解方程的一阶导数 def diff_funtion():         fh = function()     dth = fh.diff() return dth # 定义牛顿法 def Newton(x0,eps,maxiter):     # x1 = x0 - f(x0)/f`(x0)  ==>  x2 = x1 - f(x1)/f`(x1)     x = symbols('x')  # 符号变量的定义     fh = function() #方程     dfh = diff_funtion() #方程一阶导     x_n = x0 # x_next代表x(n+1)       print('%5s %14s %23s' %('迭代次数','计算结果','误差'))     #利用牛顿法逐渐逼近精确解     for k in range(maxiter):         x_b = x_n # x_before代表x(n)         fx = fh.evalf(subs = {x:x_b}) # 方程在xn处的数值         dfx = dfh.evalf(subs = {x:x_b}) # 方程的一阶导在xn处的数值         x_n = x_b - fx/dfx         errval = abs(fh.evalf(subs = {x:x_n})) #第k次迭代误差大小           # 输出迭代过程         print('%5d %25.15f %25.15f' % (k+1,x_n,errval))         # 如果误差小于给定的精度要求，则退出         if errval < eps:             break     if k+1 <= maxiter-1:         print('方程在满足精度' + str(eps) + '的条件下，近似解为：'               + str(x_n) + '，误差是：' + str(errval))     else:         print('牛顿迭代法求解数值逼近，已达到最大迭代次数，可能不收敛或精度过高...')   fh = function() plot(fh) #画出函数图像   x0 = float(input('请输入迭代初始值：')) eps = float(input('请输入方程解的精度要求：')) maxiter = int(input('请输入最大迭代次数：')) Newton(x0,eps,maxiter)

运行结果如下：

用plot()画出函数图像方便选择迭代初始点

结语

这里简单介绍了牛顿法，其本质是函数泰勒展开后，取一阶来近似，即用线性函数近似f(x)，毕竟线性函数更好算。应用牛顿法我们几乎可以求解所有方程的近似解，而且精度极高。但通过上诉求方程-e^x*sin(x)+cos(x)-1=0的解我们可以看到牛顿法不能一次性求出所有的解，即当方程有多个解时，选择不同的初始迭代点得到方程的解就不同。

作者：唐雷清

实习编辑：李欣容

稿件来源：深度学习与文旅应用实验室（DLETA）