• 背景：如何使用python求解多元多次方程组或者非线性方程组。 原创内容，转载注明出处！请勿用于商业用途！ （上篇用python拟合2019nCov感染人数的文章被不少博主转载了，发的比较早，不少博主在文章基础上添加新内容...
• 每次求解都有错误： 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 \...

%%代码：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
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(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;
> 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感染人数的文章被不少博主转载了，发的比较早，不少博主在文章基础上添加新内容也...

背景：
如何使用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!")
博主其他文章：

展开全文
• 多元非线性方程组求解(牛顿迭代法,含matlab代码)
• fslove - Matlab 求解多元多次方程组简介: 之前看到网上的一些资料良莠不齐, 各种转载之类的, 根本无法解决实际问题, 所以我打算把自己的学到的总结一下, 以实例出发讲解 fsolve.示例如下:$\begin{cases} 2x_1 - x...  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

展开全文
• 解二元一次方程组（可以扩展至解多元多次方程组）解二元一次方程组（可以扩展至解多元多次方程组）解二元一次方程组（可以扩展至解多元多次方程组）解二元一次方程组（可以扩展至解多元多次方程组
• Wegstein法注意事项 应注意如果x1和x2两点选择不当则连线的斜率等于1与直线y=x无交点从而迭代无法进行这就是Wegstein法应当避免的陷井引入一个量C Wegstein法注意事项 令q1-C 当q0时Wegstein法退化为简单的不动点...
• //经典牛顿迭代法C++实现#include#include#define N 3 // 非线性方程组中方程个数、未知量个数#define Epsilon 0.000000001 // 差向量1范数的上限#define Max 1000000000000 //最大迭代次数using namespace std;...

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

展开全文
• fslove - Matlab求解多元多次方程组简介： 之前看到网上的一些资料良莠不齐，各种转载之类的，根本无法解决实际问题，所以我打算把自己的学到的总结一下，以实例出发讲解fsolve。示例如下：{ 2 x 1 − x 2 = e a x 1...
• 求解方程： 其中，已知量为： x_TE_1=1.23; y_TE_1=3.3232; z_TE_1=0.9876; h_T=0; R_E=6378137; R_P=6356752; x_ES=1; y_ES=1; z_ES=1; 求解：x_TE，y_TE，z_TE 代码： x_TE_1=1.23;y_TE_1=3.3232; z_...
• excel求解二元一次、或者多元多次方程组，可以用规划求解，这是大约17年我了解到的方法，也在工作生活中多次用到。 但是书里的方法大不一样。线性代数里的知识在excel里大展拳脚，我当时仿佛看到了新世界。跟着步骤...
• #from sympy import * from scipy.optimize import * import numpy as np import pylab as pl #x = symbols('x') #y = symbols('y') #x_set = [100, 86, 20] #y_set = [80, 40 ,60] #w_set = [425, 320, 220] ...
• 多元一次不定方程的强力算法---同余筛数法uniqueleion6152018-10-24Python100例——第五章----不定方程的解wdt338516542013-07-19高斯消元求解多元次方程组nikelong021632016-03-26二元一次不定方程的快速解法...
• 三个未知量构成一个方程式，该CSV文件中一共有N行数据有关[x, y, z]的系数，求解三个未知量[x, y, z]的值。 文章目录 前言 一、工具包 二、使用步骤 1.读入文件 2.编写方程 总结 前言 三个未知量[x,...
• python 求解三元一次方程组，三元一次方程组为： k00+k11*2+k22*5=11 k00+k11*7+k22*6=2 k00+k11*6+k22*9=7 demo： from sympy import * k00 = Symbol('k00') k11 = Symbol('k11') k22 = Symbol('k22') # 解三元...
• 例如解：x5+x4+x3+x2+x=7.316 from sympy import * x = symbols('x') print(solve(Eq(x**5+x**4+x**3+x**2+x,7.316),x)) print(1.12968077101867**5+1.12968077101867**4+1.12968077101867**3+1.12968077101867**2+...
• 我想求解一个非常复杂的方程组，原理很简单，就是解三个方程求出三个未知数a、b、c，但每一个方程都包含exp、log、sqrt函数，不知道能否用MATLAB求出解来，我试着写了程序，但是运行不出来，我第一用MATLAB，不...
• clear allclose allsyms pi000 pi001 pi002 pi003 pi010 pi011 pi012 pi013 pi020 pi021 pi022 pi023 pi100 pi101 pi102 pi103 pi110 pi111 pi112 pi113 pi120 pi121 pi122 pi123 a1 a2 a3 b11 b21 b22 b31 b32 b33e...
• 本文列出了使用Excel中解多元次方程组的三种方法：矩阵解法、用克莱姆法则和用规划求解的方法。方法一：矩阵解法原理：对于由n个未知数，n个方程组成的多元次方程组：写成矩阵形式为Ax=b，其中A为系数n*n方阵，x...
• 本文实例讲述了C++实现的求解多元次方程。分享给大家供大家参考，具体如下：注：这里计算的是n*n的等距矩阵，代码如下：#include#include#include#includeusing namespace std;void print(double (*pArray)[4], ...
• 吾觉得要手工计算代入,求出一元高次方程,再用roots求解,或者直接用solve求解.[x1,x2,x3]=solve('x1+x2*x3+x3=100','x1-x2+2*x3=90','x1*x3+x2*x3=300')结果是：x1 =-(1/6*(558900+60*i*6884535^(1/2))^(1/3)+1160/...
• C语言解多元次方程组(矩阵法)发布时间:2016年06月08日 评论数：抢沙发阅读数：2560#include #include #define Han 200//(可自设)多元次方程组有n行n+1列(的一列是等号右边的值)，给出行数就能确定矩阵，#...
• 分析：这一道题是四元一次方程，存在两个限制条件：1是要求各种书最少买一本，2是最多剩余2元。那么我们可以先每一种书各买一本，花掉3+5+7+11=26元，还剩44元，这44元可以任意分配，这样就解决了第一个限制条件，...

...