%%代码：m文件代码：
 
function eq=myfun(Vd0,Vd1,Vd2,Vd3,Vd4,Vc2,Vc3,Vc4,Vc5,Vb1,Vb2,Vb3,Vb4,Pd1,Pd2,Pd3,Pd4,Pc1,Pc2,Pc3,Pc4)
 
syms Vd0 Vd1 Vd2 Vd3 Vd4 Vc2 Vc3 Vc4 Vc5 Vb1 Vb2 Vb3 Vb4 Pd1 Pd2 Pd3 Pd4 Pc1 Pc2 Pc3 Pc4
 
%%global Vd0 Vd1 Vd2 Vd3 Vd4 Vc2 Vc3 Vc4 Vc5 Vb1 Vb2 Vb3 Vb4 Pd1 Pd2 Pd3 Pd4 Pc1 Pc2 Pc3 Pc4
 
g=9.8
 
de=0.25
 
Cd=0.9
 
Cc=0.2
 
P0=100000
 
Ab=48
 
Ad=60
 
Ac=60
 
H=1.5
 
Kd=1.5
 
Kc=2.5
 
Kb=0.2
 
N=4
 
Vc1=0
 
Vc0=0
 
Vb0=0
 
Pd0=101000
 
Pc5=101000
 
eq(1)=Vc1+(Ab/Ac)^2*Vb1-Vc2;
 
eq(2)=Vc2+(Ab/Ac)^2*Vb2-Vc3;
 
eq(3)=Vc3+(Ab/Ac)^2*Vb3-Vc4;
 
eq(4)=Vc4+(Ab/Ac)^2*Vb4-Vc5;
 
eq(5)=Vd0-(Ab/Ad)^2*Vb0-Vd1;
 
eq(6)=Vd1-(Ab/Ad)^2*Vb1-Vd2;
 
eq(7)=Vd2-(Ab/Ad)^2*Vb2-Vd3;
 
eq(8)=Vd3-(Ab/Ad)^2*Vb3-Vd4;
 
eq(9)=(Ab/Ad)^2*Vb4-Vd4;
 
eq(10)=Vc5-Vd0;
 
eq(11)=0.5*de*Vd0^2*Kd+de*(Vd1^2-Vd0^2)*(1-0.5*Cd)-Pd0+Pd1;
 
eq(12)=0.5*de*Vd1^2*Kd+de*(Vd2^2-Vd1^2)*(1-0.5*Cd)-Pd1+Pd2;
 
eq(13)=0.5*de*Vd2^2*Kd+de*(Vd3^2-Vd2^2)*(1-0.5*Cd)-Pd2+Pd3;
 
eq(14)=0.5*de*Vd3^2*Kd+de*(Vd4^2-Vd3^2)*(1-0.5*Cd)-Pd3+Pd4;
 
eq(15)=0.5*de*Vc1^2*Kc+de*(Vc2^2-Vc1^2)*(1-0.5*Cc)-Pc1+Pc2;
 
eq(16)=0.5*de*Vc2^2*Kc+de*(Vc3^2-Vc2^2)*(1-0.5*Cc)-Pc2+Pc3;
 
eq(17)=0.5*de*Vc3^2*Kc+de*(Vc4^2-Vc3^2)*(1-0.5*Cc)-Pc3+Pc4;
 
eq(18)=0.5*de*Vc4^2*Kc+de*(Vc5^2-Vc4^2)*(1-0.5*Cc)-Pc4+Pc5;
 
eq(19)=0.5*de*Vb1^2*Kb+de*g*H-Pd1+Pc1;
 
eq(20)=0.5*de*Vb2^2*Kb+de*g*H-Pd2+Pc2;
 
eq(21)=0.5*de*Vb3^2*Kb+de*g*H-Pd3+Pc3;
 
eq(22)=0.5*de*Vb4^2*Kb+de*g*H-Pd4+Pc4;
 
end
 
%%
 
运行代码：
 
clear
 
clc
 
x0=[0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,100000,100000,100000,100000,100000,100000,100000,100000]
 
fsolve(@myfun,x0,optimset(\'Display\',\'on\'))
 
XX=[Vd0 Vd1 Vd2 Vd3 Vd4 Vc2 Vc3 Vc4 Vc5 Vb1 Vb2 Vb3 Vb4 Pd1 Pd2 Pd3 Pd4 Pc1 Pc2 Pc3 Pc4]\'
 
%%
 
一共21个未知数。每次求解都有错误：
 
Warning: Trust-region-dogleg algorithm of FSOLVE cannot handle non-square systems;
 
using Levenberg-Marquardt algorithm instead.
 
> In fsolve at 324
 
??? Undefined function or method \'full\' for input arguments of type \'sym\'.
 
Error in ==> levenbergMarquardt at 67
 
costFun = full(costFun);
 
Error in ==> fsolve at 385
 
[x,FVAL,JACOB,EXITFLAG,OUTPUT,msgData] = ...
 
%%
 
已经改了好久了每次都有各种各样的问题，希望大神能帮我解决问题到能算出结果。。。

 
背景：
 
如何使用python求解多元多次方程组或者非线性方程组。
 
原创内容，转载注明出处！请勿用于商业用途！
 
(上篇用python拟合2019nCov感染人数的文章被不少博主转载了，发的比较早，不少博主在文章基础上添加新内容也新发了新的更新后的预测，或者加入一些新的模块。博文链接如下：)
 
目录
 
一、多元多次方程
 
1.1 定义
 
我们常见的方程组有一元一次方程组，比如x+3=5这种，很简单很好解。
 
二元一次方程组，即方程组中有两个未知数，未知数的最高次数为1.
 
二元二次方程组：方程组中有两个未知数，未知数的最高次数为2.。此类方程组均有公式解法或者成形的解法。
 
但是面临多元多次方程组，解法错综复杂，是数学家们研究的内容。为了更好的解决此类问题，我们可以用python来实现。
 
1.2 例子
 
多元多次方程组例如下面这种，三元二次方程组：
 

 
下面这种，二元二次方程组。
 

 
第二个方程组实在比较复杂，因此需要借助python。
 
二、python求解工具包
 
python求解方程组的工具包较多。例如：
 
numpy：numpy.linalg.solve 可以直接求解线性方程组，numpy是python非常常用的包，解的方程也较为初级。
 
scipy：from scipy.optimize import fsolve，可以求解非线性方程组，使用较为方便，但是解集并不完备，可能漏掉一下解(后文会给个例子)scipy可以用于数学、科学、工程领域的常用软件包,可以处理插值、积分、优化，相对较初级易用
 
sympy：此工具包功能相对强大，支持符号计算、高精度计算、解方程、微积分、组合数学、离散数学、几何学、概率与统计、物理学等方面的功能。github地址：
 
https://github.com/sympy/sympy
 
sage，不支持位运算，z3约束求解器，等其他工具包，本文不详述，感兴趣的可以查找相应的内容。
 
本文详细讲述scipy以及sympy求解多次方程的方法。
 
三、scipy方法
 
3.1 使用scipy的fsolve求解
 
我们只将求解方程的部分。
 
用fsolve相对初级，也相对简单易操作，代码较为简单，只用将方程的表达式写出运行即可。fsolve近似看作用最小二乘法求解。不够很强大，很多情况下解集不完备或者无法解出。
 
例如对于
 

 
，首先要定义相应的函数：
 
def solve_function(unsolved_value):
 
x,y,z=unsolved_value[0],unsolved_value[1],unsolved_value[2]
 
return [
 
x**2+y**2-10,
 
y**2+z**2-34,
 
x**2+z**2-26,
 
]
 
求解函数三个公式都为0时候的解，中括号内为初值[0, 0, 0]
 
solved=fsolve(solve_function,[0, 0, 0])
 
全部代码：
 
#!/usr/bin/python
 
# -*- coding: UTF-8 -*-
 
"""
 
python解方程
 
"""
 
from scipy.optimize import fsolve
 
def solve_function(unsolved_value):
 
x,y,z=unsolved_value[0],unsolved_value[1],unsolved_value[2]
 
return [
 
x**2+y**2-10,
 
y**2+z**2-34,
 
x**2+z**2-26,
 
]
 
solved=fsolve(solve_function,[0, 0, 0])
 
print(solved)
 
print("Program done!")
 
"""
 
运行结果：
 
[-1. 3. 5.]
 
Program done!
 
"""
 
看出运行结果来看，此结果并非完备解集。因为x，y，z都是可正可负。例如1或者-1，3或者-3，5或者-5，但是此工具包只能解出一个解。
 
3.2 非完备解
 
显而易见，x**2-9=0的解为3或者-3
 
def solve_function(unsolved_value):
 
x=unsolved_value[0]
 
return [
 
x**2-9,
 
]
 
solved=fsolve(solve_function,[0])
 
但是程序只能得出一个结果3，但是得不到-3
 
3.3 非线性方程的解
 
最简单的sin(x)=0.5,则x可能为π/6或者 5π/6
 
#!/usr/bin/python
 
# -*- coding: UTF-8 -*-
 
"""
 
python解方程
 
"""
 
from scipy.optimize import fsolve
 
from math import sin,cos
 
def solve_function(unsolved_value):
 
x=unsolved_value[0]
 
return [
 
sin(x)-0.5
 
]
 
solved=fsolve(solve_function,[3.14])
 
print(solved)
 
solved=fsolve(solve_function,[0])
 
print(solved)
 
print("Program done!")
 
运行结果为：
 
[2.61799388]
 
[0.52359878]
 
Program done!
 
可以解出π/6或者 5π/6，中括号内为初始迭代的值。
 
3.4 无法求解
 
部分较难情况无法求解
 

 
#!/usr/bin/python
 
# -*- coding: UTF-8 -*-
 
"""
 
python解方程
 
"""
 
from scipy.optimize import fsolve
 
def solve_function(unsolved_value):
 
x,y=unsolved_value[0],unsolved_value[1]
 
return [
 
x*x+2*x*y,
 
2*x*y-2*y*y
 
]
 
solved=fsolve(solve_function,[6, -3])
 
print(solved)
 
print("Program done!")
 
无法求解会给出报错，和用最小二乘法迭代得到明显错误的解。
 
[1.64526700e-115 1.33665018e-115]
 
A:\python\python\lib\site-packages\scipy\optimize\minpack.py:162: RuntimeWarning: The number of calls to function has reached maxfev = 600.
 
Program done!
 
warnings.warn(msg, RuntimeWarning)
 
四、sympy工具包求解
 
没安装可以在teiminal中pip install sympy，此工具包涉及支持符号计算、高精度计算、模式匹配、绘图、解方程、微积分、组合数学、离散 数学、几何学、概率与统计、物理学等方面的功能。功能较为强大，解方程组时性能也较好。
 
4.1 二元一次方程组
 

 
较为简单，
 
from sympy import *
 
# 二元一次方程
 
x = Symbol('x')
 
y = Symbol('y')
 
solved_value=solve([2*x+y-1, x-2*y], [x, y])
 
print(solved_value)
 
此方法较为简单，但是相应的自变量应当写成符号的形式，x=Symbol('x')
 
求解后有分数解：
 
{x: 2/5, y: 1/5}
 
Program done!
 
4.2 多解
 
多解情况与复数解
 
例如，多个解的情况，sympy可以很好的进行求解
 

 
x = Symbol('x')
 
solved_value=solve([x**2-9], [x])
 
print(solved_value)
 
输出结果：
 
[(-3,), (3,)]
 
4.3 复数解
 
复数解也可以很好解出：
 

 
# 复数解
 
solved_value = solve([x ** 2 + 9], [x])
 
print(solved_value)
 
solved_value = solve([x ** 4 - 9], [x])
 
print(solved_value)
 
"""
 
运行结果：
 
[(-3*I,), (3*I,)]
 
[(-sqrt(3),), (sqrt(3),), (-sqrt(3)*I,), (sqrt(3)*I,)]
 
"""
 
复数解也能较好解出
 
4.4 非线性求解
 
比如三角函数：
 

 
程序均能较好解出
 
# 非线性解
 
solved_value = solve([sin(x) - 0.5], [x])
 
print(solved_value)
 
solved_value = solve([sin(x) - 1], [x])
 
print(solved_value)
 
"""
 
[(0.523598775598299,), (2.61799387799149,)]
 
[(pi/2,)]
 
"""
 
4.5 较为复杂的二元二次方程
 

 
此题较难，无论人来算，很难算出，用scipy工具包也迭代不出解。但是sympy强大的功能可以很好的解出此方程。
 
# 二元二次方程组
 
x = Symbol('x')
 
y= Symbol('y')
 
solved_value=solve([x**2+2*x*y-6,2*x*y-2*y**2+3], [x,y])
 
print(solved_value)
 
有四组实数解：
 
[(-(-3 + sqrt(13))*sqrt(sqrt(13)/2 + 2), -sqrt(sqrt(13)/2 + 2)),
 
((-3 + sqrt(13))*sqrt(sqrt(13)/2 + 2), sqrt(sqrt(13)/2 + 2)),
 
(-sqrt(2 - sqrt(13)/2)*(-sqrt(13) - 3), -sqrt(2 - sqrt(13)/2)),
 
(sqrt(2 - sqrt(13)/2)*(-sqrt(13) - 3), sqrt(2 - sqrt(13)/2))]
 

 
复杂的问题终于解出，有四组实数解！
 
五、全部代码
 
#!/usr/bin/python
 
# -*- coding: UTF-8 -*-
 
"""
 
python解方程
 
created by xingxinagrui on 2020.2.24
 
"""
 
from scipy.optimize import fsolve
 
from math import sin,cos
 
from sympy import *
 
# 1-4 scipy
 
# 5-7 sympy
 
part=7
 
if part==1:
 
# 求解非线性方程组
 
def solve_function(unsolved_value):
 
x=unsolved_value[0]
 
return [
 
sin(x)-0.5
 
]
 
solved=fsolve(solve_function,[3.14])
 
print(solved)
 
solved=fsolve(solve_function,[0])
 
print(solved)
 
if part==2:
 
# 求解三元二次方程组
 
def solve_function(unsolved_value):
 
x, y, z = unsolved_value[0], unsolved_value[1], unsolved_value[2]
 
return [
 
x ** 2 + y ** 2 - 10,
 
y ** 2 + z ** 2 - 34,
 
x ** 2 + z ** 2 - 26,
 
]
 
solved = fsolve(solve_function, [0, 0, 0])
 
print(solved)
 
if part==3:
 
#解的非完备性
 
def solve_function(unsolved_value):
 
x = unsolved_value[0]
 
return [
 
x ** 2 - 9,
 
]
 
solved = fsolve(solve_function, [0])
 
print(solved)
 
if part == 4:
 
# 较难无法求解
 
def solve_function(unsolved_value):
 
x, y = unsolved_value[0], unsolved_value[1]
 
return [
 
x * x + 2 * x * y,
 
2 * x * y - 2 * y * y
 
]
 
solved = fsolve(solve_function, [6, -3])
 
print(solved)
 
if part == 5:
 
# 二元一次方程
 
x = Symbol('x')
 
y = Symbol('y')
 
solved_value=solve([2*x+y-1, x-2*y], [x, y])
 
print(solved_value)
 
if part == 6:
 
# 多解情况
 
x = Symbol('x')
 
solved_value=solve([x**2-9], [x])
 
print(solved_value)
 
# 复数解
 
solved_value = solve([x ** 2 + 9], [x])
 
print(solved_value)
 
solved_value = solve([x ** 4 - 9], [x])
 
print(solved_value)
 
# 非线性解
 
solved_value = solve([sin(x) - 0.5], [x])
 
print(solved_value)
 
solved_value = solve([sin(x) - 1], [x])
 
print(solved_value)
 
if part == 7:
 
# 二元二次方程组
 
x = Symbol('x')
 
y= Symbol('y')
 
solved_value=solve([x**2+2*x*y-6,2*x*y-2*y**2+3], [x,y])
 
print(solved_value)
 
print("Program done!")
 
博主其他文章：

 
fslove - Matlab 求解多元多次方程组
 
简介: 之前看到网上的一些资料良莠不齐, 各种转载之类的, 根本无法解决实际问题, 所以我打算把自己的学到的总结一下, 以实例出发讲解 fsolve.
 
示例如下:
 
\[ \begin{cases} 2x_1 - x_2 = e^{ax_1} \-x_1 + 2x_2 = e^{ax_2} \\end{cases} \]
 
具体的求解过程在后面 点击跳转
 
1. fsolve 的基本使用
 
调用格式一:
 
X = fslove(FUN,X0)
 
功能: 给定初值 X0, 求解方程组的解, X 就是返回的解
 
调用格式二:
 
X = fsolve(FUN,X0,OPTIONS)
 
功能: 同上, 并解决默认参数优化为 options 指定值
 
调用格式三:
 
[X,FVAL] = fslove(FUN,X0,...)
 
功能: 返回 X 处目标函数值
 
调用格式四:
 
[X,FVAL,EXITFLAG] = fslove(FUN,X0,...)
 
功能: 返回 EXITFLAG 的值, 用来描述计算退出的条件, 其中 EXITFLAG 取值和相应的含义如下表.(主要作为判断条件来使用)EXITFLAG含义1函数 fslove 收敛于解 X 处
 
2X 的变化小于限制
 
3残差变化小于限制
 
4重要搜索方向小于限制
 
0达到最大迭代次数或者评价标准
 
-1算法由输出函数终止
 
-2算法无法收敛到解的点
 
-3信赖域半径太小
 
-4线搜索在当前不能充分减少残差
 
调用格式五:
 
[X,FVAL,EXITFLAG,OUTPUT] = fslove(FUN,X0,...)
 
功能: 包含 OUTPUT 的输出
 
调用格式六:
 
[X,FVAL,EXITFLAG,OUTPUT,JACOB] = fslove(FUN,X0,...)
 
功能: 返回雅各比矩阵
 
2. 方程求解
 
(1) 编制函数文件 fun.m
 
编写函数主要用来书写函数的表达式.
 
function f = fun(x,a,b,c) % b c 可以是随意的参数f1=2*x(1)-x(2)-exp(a*x(1));
 
f2=-x(1)+2*x(2)-exp(a*x(2));
 
f=[f1;f2];
 
% 也可以写成下面的方式
 
% f = [2*x(1)-x(2)-exp(a*x(1));-x(1)+2*x(2)-exp(a*x(2))];
 
(2) 给定函数的参数值和初值 (解在周围寻找)
 
调用求解函数 fslove>>a=-1;
 
>>x0=[-5,-4];
 
>>[x,FVAL,EXITFLAG,OUTPUT,JACOB]=fsolve(@(x)fun(x,a,1,1),x0);
 
@(x)fun(x,a,1,1) 调用 fun 函数, 函数的参数是 a,1,1, 求解 x 的值
 
执行后调用 x 返回, 也就是 X 的解.x=
 
??0.5671???0.5671
 
调用 FVAL 显示在目标解的函数值, 可以看出, FVAL 越小越接近真实解.FVAL=
 
??1.0e-09*
 
??-0.4242
 
??-0.3753
 
调用 EXITFLAG 结合上面的表格可以知道, 函数 FSOLVE 收敛于解 X 处.
 
EXITFLAG =
 
??1
 
来源: http://www.bubuko.com/infodetail-3213047.html

 
//经典牛顿迭代法C++实现
 
#include
 
#include
 
#define N 3 // 非线性方程组中方程个数、未知量个数
 
#define Epsilon 0.000000001 // 差向量1范数的上限
 
#define Max 1000000000000 //最大迭代次数
 
using namespace std;
 
const int N2=2*N;
 
int main()
 
{
 
void ff(double xx[N],double yy[N]); //计算向量函数的因变量向量yy[N]
 
void ffjacobian(double xx[N],double yy[N][N]); //计算雅克比矩阵yy[N][N]
 
void inv_jacobian(double yy[N][N],double inv[N][N]); //计算雅克比矩阵的逆矩阵inv
 
void newdundiedai(double x0[N], double inv[N][N],double y0[N],double x1[N]); //由近似解向量 x0 计算近似解向量 x1
 
// double x0[N]={82.7995,-5.92913,361.667};
 
double x0[N]={-5,-6,-7};
 
double y0[N],jacobian[N][N],invjacobian[N][N],x1[N],errornorm;
 
int i,j,iter=0;
 
cout<
 
for (i=0;i
 
cout<
 
cout<
 
do
 
{
 
iter=iter+1;
 
cout<
 
//计算向量函数的因变量向量 y0
 
ff(x0,y0);//double aaa = exp(-5);
 
//计算雅克比矩阵 jacobian
 
ffjacobian(x0,jacobian);
 
//计算雅克比矩阵的逆矩阵 invjacobian
 
inv_jacobian(jacobian,invjacobian);
 
//由近似解向量 x0 计算近似解向量 x1
 
newdundiedai(x0, invjacobian,y0,x1);
 
//计算差向量的1范数errornorm
 
errornorm=0;
 
for (i=0;i
 
errornorm=errornorm+fabs(x1[i]-x0[i]);
 
if (errornorm
 
for (i=0;i
 
x0[i]=x1[i];
 
} while (iter
 
return 0;
 
}
 
void ff(double xx[N],double yy[N])
 
{
 
double x,y,z;
 
int i;
 
x=xx[0];
 
y=xx[1];
 
z=xx[2];
 
// yy[0]=3*x-2*y+4*z-11;
 
// yy[1]=2*x*x+2*x+y*y+3*y+4*z*z-27;
 
// yy[2]=x+2*y+3*z-14;
 
yy[0]=3*x-cos(y*z)-0.5;
 
yy[1]=x*x-81*(y+0.1)*(y+0.1)+sin(z)+1.06;
 
yy[2]=exp(-x*y)+20*z+10.0/3.0*3.14159-1;
 
cout<
 
for( i=0;i
 
cout<
 
cout<
 
cout<
 
}
 
void ffjacobian(double xx[N],double yy[N][N])
 
{
 
double x,y,z;
 
int i,j;
 
x=xx[0];
 
y=xx[1];
 
z=xx[2];
 
//jacobian have n*n element
 
// yy[0][0]=2*x-2;
 
// yy[0][1]=-1;
 
// yy[1][0]=2*x;
 
// yy[1][1]=8*y;
 
// yy[0][0]=3;
 
// yy[0][1]=-2;
 
// yy[0][2]=4;
 
// yy[1][0]=4*x+2;
 
// yy[1][1]=2*y+3;
 
// yy[1][2]=8*z;
 
// yy[2][0]=1;
 
// yy[2][1]=2;
 
// yy[2][2]=3;
 
yy[0][0]=3;
 
yy[0][1]=sin(y*z)*z;
 
yy[0][2]=0;
 
yy[1][0]=2*x;
 
yy[1][1]=-81*2*(y+0.1);
 
yy[1][2]=cos(z);
 
yy[2][0]=exp(-x*y)*(-y);
 
yy[2][1]=exp(-x*y)*(-x);
 
yy[2][2]=20;
 
cout<
 
for( i=0;i
 
{for(j=0;j
 
cout<
 
cout<
 
}
 
cout<
 
}
 
void inv_jacobian(double yy[N][N],double inv[N][N])
 
{
 
double aug[N][N2],L;
 
int i,j,k;
 
cout<
 
for (i=0;i
 
{ for(j=0;j
 
aug[i][j]=yy[i][j];
 
for(j=N;j
 
if(j==i+N) aug[i][j]=1;
 
else aug[i][j]=0;
 
}
 
for (i=0;i
 
{ for(j=0;j
 
cout<
 
cout<
 
}
 
cout<
 
for (i=0;i
 
{
 
for (k=i+1;k
 
{L=-aug[k][i]/aug[i][i];
 
for(j=i;j
 
aug[k][j]=aug[k][j]+L*aug[i][j];
 
}
 
}
 
for (i=0;i
 
{ for(j=0;j
 
cout<
 
cout<
 
}
 
cout<
 
for (i=N-1;i>0;i--) -0
 
{
 
for (k=i-1;k>=0;k--)
 
{L=-aug[k][i]/aug[i][i];
 
for(j=N2-1;j>=0;j--)
 
aug[k][j]=aug[k][j]+L*aug[i][j];
 
}
 
}
 
for (i=0;i
 
{ for(j=0;j
 
cout<
 
cout<
 
}
 
cout<
 
for (i=N-1;i>=0;i--)
 
for(j=N2-1;j>=0;j--)
 
aug[i][j]=aug[i][j]/aug[i][i];
 
for (i=0;i
 
{ for(j=0;j
 
cout<
 
cout<
 
for(j=N;j
 
inv[i][j-N]=aug[i][j];
 
}
 
cout<
 
cout<
 
for (i=0;i
 
{ for(j=0;j
 
cout<
 
cout<
 
}
 
cout<
 
}
 
void newdundiedai(double x0[N], double inv[N][N],double y0[N],double x1[N])
 
{
 
int i,j;
 
double sum=0;
 
for(i=0;i
 
{ sum=0;
 
for(j=0;j
 
sum=sum+inv[i][j]*y0[j];
 
x1[i]=x0[i]-sum;
 
}
 
cout<
 
for (i=0;i
 
cout<
 
cout<
 
}
 
牛顿法的算法步骤：
 

 

 

 
其实看懂了就很简单：就是要求两个矩阵。雅克比矩阵及其逆矩阵。雅克比矩阵就是上图，雅克比矩阵就是偏导数，求完雅克比矩阵，矩阵里面元素都是常数了，就是再去求这个常数矩阵的逆矩阵，求逆矩阵代码一大把。
 
再给出刚开始代码链接里面一个二元二次方程组的例子的源代码：
 
//经典牛顿迭代法C++实现
 
#include
 
#include
 
#define N 2 // 非线性方程组中方程个数、未知量个数
 
#define Epsilon 0.0001 // 差向量1范数的上限
 
#define Max 100 //最大迭代次数
 
using namespace std;
 
const int N2=2*N;
 
int main()
 
{
 
void ff(float xx[N],float yy[N]); //计算向量函数的因变量向量yy[N]
 
void ffjacobian(float xx[N],float yy[N][N]); //计算雅克比矩阵yy[N][N]
 
void inv_jacobian(float yy[N][N],float inv[N][N]); //计算雅克比矩阵的逆矩阵inv
 
void newdundiedai(float x0[N], float inv[N][N],float y0[N],float x1[N]); //由近似解向量 x0 计算近似解向量 x1
 
float x0[N]={12.0,10.25},y0[N],jacobian[N][N],invjacobian[N][N],x1[N],errornorm;
 
int i,j,iter=0;
 
//如果取消对x0的初始化，撤销下面两行的注释符，就可以由键盘向x0读入初始近似解向量
 
//for( i=0;i
 
// cin>>x0[i];
 
cout<
 
for (i=0;i
 
cout<
 
cout<
 
do
 
{
 
iter=iter+1;
 
cout<
 
//计算向量函数的因变量向量 y0
 
ff(x0,y0);
 
//计算雅克比矩阵 jacobian
 
ffjacobian(x0,jacobian);
 
//计算雅克比矩阵的逆矩阵 invjacobian
 
inv_jacobian(jacobian,invjacobian);
 
//由近似解向量 x0 计算近似解向量 x1
 
newdundiedai(x0, invjacobian,y0,x1);
 
//计算差向量的1范数errornorm
 
errornorm=0;
 
for (i=0;i
 
errornorm=errornorm+fabs(x1[i]-x0[i]);
 
if (errornorm
 
for (i=0;i
 
x0[i]=x1[i];
 
} while (iter
 
return 0;
 
}
 
void ff(float xx[N],float yy[N])
 
{float x,y;
 
int i;
 
x=xx[0];
 
y=xx[1];
 
yy[0]=x*x-2*x-y+0.5;
 
yy[1]=x*x+4*y*y-4;
 
cout<
 
for( i=0;i
 
cout<
 
cout<
 
cout<
 
}
 
void ffjacobian(float xx[N],float yy[N][N])
 
{
 
float x,y;
 
int i,j;
 
x=xx[0];
 
y=xx[1];
 
//jacobian have n*n element
 
yy[0][0]=2*x-2;
 
yy[0][1]=-1;
 
yy[1][0]=2*x;
 
yy[1][1]=8*y;
 
cout<
 
for( i=0;i
 
{for(j=0;j
 
cout<
 
cout<
 
}
 
cout<
 
}
 
void inv_jacobian(float yy[N][N],float inv[N][N])
 
{float aug[N][N2],L;
 
int i,j,k;
 
cout<
 
for (i=0;i
 
{ for(j=0;j
 
aug[i][j]=yy[i][j];
 
for(j=N;j
 
if(j==i+N) aug[i][j]=1;
 
else aug[i][j]=0;
 
}
 
for (i=0;i
 
{ for(j=0;j
 
cout<
 
cout<
 
}
 
cout<
 
for (i=0;i
 
{
 
for (k=i+1;k
 
{L=-aug[k][i]/aug[i][i];
 
for(j=i;j
 
aug[k][j]=aug[k][j]+L*aug[i][j];
 
}
 
}
 
for (i=0;i
 
{ for(j=0;j
 
cout<
 
cout<
 
}
 
cout<
 
for (i=N-1;i>0;i--)
 
{
 
for (k=i-1;k>=0;k--)
 
{L=-aug[k][i]/aug[i][i];
 
for(j=N2-1;j>=0;j--)
 
aug[k][j]=aug[k][j]+L*aug[i][j];
 
}
 
}
 
for (i=0;i
 
{ for(j=0;j
 
cout<
 
cout<
 
}
 
cout<
 
for (i=N-1;i>=0;i--)
 
for(j=N2-1;j>=0;j--)
 
aug[i][j]=aug[i][j]/aug[i][i];
 
for (i=0;i
 
{ for(j=0;j
 
cout<
 
cout<
 
for(j=N;j
 
inv[i][j-N]=aug[i][j];
 
}
 
cout<
 
cout<
 
for (i=0;i
 
{ for(j=0;j
 
cout<
 
cout<
 
}
 
cout<
 
}
 
void newdundiedai(float x0[N], float inv[N][N],float y0[N],float x1[N])
 
{
 
int i,j;
 
float sum=0;
 
for(i=0;i
 
{ sum=0;
 
for(j=0;j
 
sum=sum+inv[i][j]*y0[j];
 
x1[i]=x0[i]-sum;
 
}
 
cout<
 
for (i=0;i
 
cout<
 
cout<
 
}
 
不过，弄好了是简单，这段代码也消耗了我一天的时间，首先是看原理，因为上学的时候只记得电力系统计算潮流那段有个叫雅克比矩阵的东东的，至于是干嘛用的，完全忘记了，当时只应付考试了。看完原理然后就是代码，有的代码是错的，调试半天也不行。
 
原文：https://www.cnblogs.com/yanghailin/p/11930494.html