修正牛顿拉普森法(MNR)解非线性方程组(python,数值积分)

第二十三篇 非线性方程组的修正牛顿拉普森法

在整个牛顿-拉夫森算法中，由于需要不停地求矩阵的逆，所以如果n变大，意味着每次迭代都需要大量的计算。此外，由于系数矩阵通常是非对称的，而且每次迭代都会改变系数矩阵的值，所以因式分解法也没有任何方便可言。修正的Newton-Raphson方法目的是减少每次迭代执行的计算量，但缺点是需要更多的迭代次数来实现收敛。不需要每次迭代都更新和求逆雅可比矩阵，可以周期性地更新它，比如m次迭代更新一下。
一个相对极端但简便的方法是，只计算并求解一次雅可比矩阵的逆，对应于猜测的初始值。这与之前提到过的对单根的求解策略相类似。正常情况，它需要的迭代次数比完整的牛顿-拉夫森过程更多。在实践中，迭代次数和计算复杂二者的平衡还需要进一步衡量。
计算实例：
使用修正牛顿拉普森方法去发现方程在x=-2.4和y=-0.1附近的根

首先根据下表把方程重新排列成标准形式

以上篇提到的雅可比矩阵对上面公式求微分

对小型非线性方程组，很容易取得雅可比矩阵的逆，为

当初始猜测值x1=-2.4和x2=-0.1时，

在修正牛顿拉普森方法中，[J]矩阵的逆取恒定值。
第一次迭代，初始猜测值x1=-2.4和x2=-0.1，因此

因此

第二次迭代：
更新值为

之后

因此

第三次迭代：
更新值

程序如下:
其中有一个主程序和一个检查是否收敛的子程序checkit，还有两个函数程序，分别为函数的表达式，和函数导数的表达式
主程序：

#非线性方程组的修正牛顿拉普森法
import numpy as np
import B
n=2;tol=1.0e-5;limit=100
x1=np.zeros((n,1))
inv=np.zeros((n,1))
x0=np.ones((2,1),dtype=np.float)
print('开始的猜测值',x0[:,0])
print('前几次迭代值')
iters=0
def f38(x):
    f38=np.zeros((2,1),dtype=np.float)
    f38[0,0]=2.0*x[0,0]**2+x[1,0]**2-4.32
    f38[1,0]=x[0,0]**2-x[1,0]**2
    return f38
def f38dash(x):
    f38dash=np.zeros((2,2),dtype=np.float)
    f38dash[0,0]=4.0*x[0,0]
    f38dash[1,0]=2.0*x[0,0]
    f38dash[0,1]=2.0*x[1,0]
    f38dash[1,1]=-2.0*x[1,0]
    return f38dash
inv=np.linalg.inv(f38dash(x0))
while(True):
    iters=iters+1
    x1[:,0]=x0[:,0]-np.dot(np.transpose(f38(x0)),inv)
    if  B.checkit(x1,x0,tol)==True or iters==limit:
        break
    x0[:,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.0
  converged=True
  for 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=False
  checkit=converged
  return  checkit

终端输出结果：
)
程序结果与上篇牛顿拉普森的计算实例结果一致。