第二十二篇 非线性方程组的牛顿拉普森解
这种方法是基于几个变量的泰勒展开式来求解的。例如，对于两个方程，假设我们展开一个关于根(x1,x2)的猜想的泰勒级数。对于x1方向上的“一小步”Δx1和x2方向上的“一小步”Δx2，一阶泰勒展开为
其中
假设(x1+Δx1,x2+Δx2)是方程的根，那么对应f1和f2的值就会为0，上式将会变成
或者以矩阵形式，左边为雅可比矩阵，其行列式称为“雅可比”。显然，如果雅可比为零，则迭代过程失败，如果它接近零，则收敛速度可能会较慢。
其中所有函数和导数都在初始猜测值(x1,x2)处取值，因此，对于n个非线性方程组，我们必须求解n个线性方程来求得向量的变化值
使用上面的值去更新所有的向量
这个过程不断地重复，直到上面方程的解向量基本不再改变
计算实例:
使用牛顿拉普森法去发现方程在x=1.8和y=0.5附近的根
首先我们根据下标符号将方程排列成标准形式
通过上面公式提到的微分形成雅可比矩阵
在小型的非线性方程组中，很容易得到雅可比矩阵的逆，因此
第一次迭代
初始猜测值，x1=1.8,x2=0.5,因此
因此
第二次迭代，更新值
因此
因此
第三次迭代
更新值
随着不断迭代Δx1、Δx2→0，就越来越接近根。
程序如下
分为一个主程序和一个检查多个向量收敛的子程序checkit，两个函数程序，分别为函数形式和函数的导数形式。程序中矩阵程法使用，np.dot；矩阵取逆采用np.linalg.inv
主程序：
#非线性方程组的牛顿拉普森法
import numpy as np
import B
n=2;tol=1.0e-5;limit=100
x1=np.zeros((n,1))
x0=np.ones((2,1),dtype=np.float)
print('开始的猜测值',x0[:,0])
print('前几次迭代值')
iters=0
def f37(x):f37=np.zeros((2,1),dtype=np.float)f37[0,0]=2.0*x[0,0]**2+x[1,0]**2-4.32f37[1,0]=x[0,0]**2-x[1,0]**2return f37
def f37dash(x):f37dash=np.zeros((2,2),dtype=np.float)f37dash[0,0]=4.0*x[0,0]f37dash[1,0]=2.0*x[0,0]f37dash[0,1]=2.0*x[1,0]f37dash[1,1]=-2.0*x[1,0]return f37dash
while(True):iters=iters+1x1[:,0]=x0[:,0]-np.dot(np.transpose(f37(x0)),np.linalg.inv(f37dash(x0)))if B.checkit(x1,x0,tol)==True or iters==limit:breakx0[:,0]=x1[:,0]if iters<5:print(x0[:,0])
print('迭代到收敛的次数',iters)
print('解',x0[:,0])
checkit
def checkit(loads,oldlds,tol):
#检查多个未知数的收敛neq=loads.shape[0]big=0.0converged=Truefor i in range(1,neq+1):if abs(loads[i-1,0])>big:big=abs(loads[i-1,0])for i in range(1,neq+1):if abs(loads[i-1,0]-oldlds[i-1,0])/big>tol:converged=Falsecheckit=convergedreturn checkit
终端输出结果如下
程序结果与计算结果一致