1 灵感来源
看书，偶然看到“解联立方程组”的内容，描述了如何用数组公式minverse和mmult求解二元一次方程组，惊到了。

二元一次方程组，通过消元的方式，可以手动求解，这是小学就练习了无数次的知识点。

excel求解二元一次、或者多元多次方程组，可以用规划求解，这是大约17年我了解到的方法，也在工作生活中多次用到。

但是书里的方法大不一样。线性代数里的知识在excel里大展拳脚，我当时仿佛看到了新世界。跟着步骤练手，成了。

脑子里产生了一个新的念头：解二元、三元、N元，变量个数不同，求解的excel就有多少个（单元格区域大小不一样不能通用）。怎么才能把不同变量数的多元一次方程组做到一个模板里呢？

经过摸索练习，最后终于做成了想要的样子（2020/6/1）。
2 思路分析
求解的部分已经解决，利用书里介绍的方法：

整理方程式，将带系数的变量放在等号左边，常量放在等号右边；
使用minverse函数计算系数矩阵的逆矩阵；
使用mmult函数计算系数矩阵的逆矩阵乘以常量矩阵，即得变量的解。

注：系数逆矩阵城西常量矩阵为什么就能得到变量的解，这一点我也没弄明白具体过程，大学时候学的线性代数已经完全还给老师了对不起。

关键的问题在于，如何在变量数量不等的条件下，实现以上过程。

即：如何实现大小变化的矩阵的范围的引用。（假设行/列数为1-50。）

当时尝试了不少方法（原谅我是个摸索中前进的小白，那时候对offset不熟），最后，通过offset函数实现了目的。

注：印象中也尝试了其他函数实现不同范围的引用，但是范围之外的单元格由于计算不到，会显示错误。offset没有这个问题，故最后敲定offset。
3 实现步骤

3.1 变量数单元格设置
由上分析，如果要做成适用于多元一次方程组的求解计算器，变量的个数是一个重要的参数。故首先设置了一个单元格，用来填入变量的个数，即多元一次方程组的元数。
3.2 系数、常量单元格设置

制作50×50的系数表；
将方程式的变量的系数填入；
将常量填写在系数单元格右侧（其他地方也可以，看个人喜好）。

注：此处设计的50×50系数矩阵，即解五十元一次方程组。理论上N元都可以，但是好像某个函数有参数限制，所以此处设置的50。具体哪个函数，找了一下，没找到，有兴趣的亲可以再求证了。
3.3 系数逆矩阵计算


使用minverse+offset函数，minverse计算系数逆矩阵，offset引用变量单元格的数值指定引用范围。

注：数组公式需按ctrl+shift+enter一起输入。
3.4 变量求解


使用mmult+offset函数计算系数矩阵的逆矩阵乘以常量矩阵，即得多元一次方程组变量的解。

注：数组公式需按ctrl+shift+enter一起输入。
4 总结&感想
总体来说，是线性代数+数组函数+offset一起协作解决了一个有点有趣又似乎没什么实际应用价值的问题。

不知道用到它的会是什么样的朋友，但是我相信它还是会遇到一直找它的小伙伴的~
【原创内容，引用请注明出处：https://blog.csdn.net/truly121314/article/details/107854989】

背景：
如何使用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处


2
X的变化小于限制


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

//经典牛顿迭代法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