牛顿迭代法求解二元非线性方程组,C++代码实现

整体迭代公式就是

X^{k+1}=X^{k}-J^{-1}\bullet F

上式中

X=\begin{bmatrix} x1\\ x2\end{bmatrix}F=\begin{bmatrix} f1\\ f2\end{bmatrix}

JF的雅克比矩阵,J= \begin{bmatrix} df1dx1 & df1dx2\\ df2dx1 & df2dx2 \end{bmatrix}

J^{-1}为雅克比矩阵的逆矩阵

实例:

求解\left\{\begin{matrix} x^{2}-2x-y=-0.5\\ x^{2}+4y^{2}=4\end{matrix}\right.的x和y的解。

上面求根问题可转化为F= \begin{bmatrix} x^{2}-2x-y+0.5\\ x^{2}+4y^{2}-4\end{bmatrix}=0的问题,即可用牛顿迭代法求解此二元非线性方程。

具体求解过程代码如下所示

//线性方程组中方程个数、未知量个数
#include <iostream> 
#include <cmath>
#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]={2.0,0.25},y0[N],jacobian[N][N],invjacobian[N][N],x1[N],errornorm;
	
	int i,j,iter=0;
	
	//如果取消对x0的初始化,撤销下面两行的注释符,就可以由键盘向x0读入初始近似解向量
	/*
	for( i=0;i<N;i++)
		cin>>x0[i];
	*/
	cout<<"初始近似解向量:"<<endl;
	
	for (i=0;i<N;i++) 
		cout<<x0[i]<<"	";
	
	cout<<endl<<endl;
	
	//牛顿迭代法迭代求解过程 
	do
	{
		iter=iter+1;
		for(int xx=0;xx<=30;xx++)  
			cout<<"==";
		cout<<endl;
		cout<<"第"<<iter<<"次迭代开始"<<endl; 
		
		ff(x0,y0); 
		
		//计算雅克比矩阵jacobian
		ffjacobian(x0,jacobian); 
		
		//计算雅克比矩阵的逆矩阵invjacobian
		inv_jacobian(jacobian,invjacobian); 
		
		//由近似解向量x0 计算近似解向量x1
		newdundiedai(x0, invjacobian,y0,x1); 
		
		//计算差向量的1范数errornorm errornorm=0;
		errornorm = 0;
		for (i=0;i<N;i++)
			errornorm=errornorm+fabs(x1[i]-x0[i]);
		
		if (errornorm<Epsilon) break;
		
		for (i=0;i<N;i++)
			x0[i]=x1[i];
	
	}while (iter<Max);
	
	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<<"向量函数的因变量向量是:"<<endl;
	for( i=0;i<N;i++)
		cout<<yy[i]<<"	";
	cout<<endl<<endl;
}

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<<"雅克比矩阵是:"<<endl;
	
	for( i=0;i<N;i++)
	{
		for(j=0;j<N;j++)
			cout<<yy[i][j]<<"	";
		cout<<endl;
	}
	
	cout<<endl;

}

void inv_jacobian(float yy[N][N],float inv[N][N]) 
{
	float aug[N][N2],L;
	int i,j,k;

	cout<<"开始计算雅克比矩阵的逆矩阵:"<<endl;

	for(i=0;i<N;i++)
	{
		for(j=0;j<N;j++)
			aug[i][j]=yy[i][j];
	
		for(j=N;j<N2;j++)
			if(j==i+N) aug[i][j]=1;
			else aug[i][j]=0;
	}
	
	for (i=0;i<N;i++)
	{ 
		for(j=0;j<N2;j++)
			cout<<aug[i][j]<<" ";
		cout<<endl;
	}
	
	cout<<endl;
	for (i=0;i<N;i++)
	{
		for (k=i+1;k<N;k++)
		{
			L=-aug[k][i]/aug[i][i];
			for(j=i;j<N2;j++)
				aug[k][j]=aug[k][j]+L*aug[i][j];
		}
	}
	
	for (i=0;i<N;i++)
	{
		for(j=0;j<N2;j++)
			cout<<aug[i][j]<<"	";
		cout<<endl;
	
	}
	cout<<endl;
	
	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<N;i++)
	{
		for(j=0;j<N2;j++)
			cout<<aug[i][j]<<"	";
		cout<<endl;
	}
	
	cout<<endl;
	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<N;i++)
	{ 
		for(j=0;j<N2;j++)
			cout<<aug[i][j]<<"	";
		
		cout<<endl;
		
		for(j=N;j<N2;j++)
			inv[i][j-N]=aug[i][j];
	}
	
	cout<<endl;
	
	cout<<"雅克比矩阵的逆矩阵:"<<endl;
	for (i=0;i<N;i++)
	{ 
		for(j=0;j<N;j++)
			cout<<inv[i][j]<<"	";
		cout<<endl;
	}
	
	cout<<endl;

}

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<N;i++)
	{ 
		sum=0;
		for(j=0;j<N;j++)
			sum=sum+inv[i][j]*y0[j];
		x1[i]=x0[i]-sum;
	}
	cout<<"近似解向量:"<<endl;
	for (i=0;i<N;i++)
		cout<<x1[i]<<"	";
	cout<<endl<<endl;

}

  • 3
    点赞
  • 37
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
牛顿迭代法可以用于求解二元线性方程组。假设我们要求解如下方程组: f1(x1, x2) = 0 f2(x1, x2) = 0 其中,x1和x2是未知量,f1和f2是已知的线性函数。 牛顿迭代法的思路是,从一个初始点(x1^0, x2^0)开始,通过不断迭代,使得每一步迭代后的点(x1^k, x2^k)都更加接近方程组的解。具体迭代公式如下: x1^(k+1) = x1^k - [J^-1(x1^k, x2^k) * F(x1^k, x2^k)]1 x2^(k+1) = x2^k - [J^-1(x1^k, x2^k) * F(x1^k, x2^k)]2 其中,J是雅可比矩阵,F是线性方程组的函数向量,[J^-1(x1^k, x2^k) * F(x1^k, x2^k)]表示矩阵J的逆与向量F的乘积。 具体的C++代码如下: ```c++ #include <iostream> #include<cmath> using namespace std; double f1(double x1, double x2) { return pow(x1, 2) + pow(x2, 2) - 4; } double f2(double x1, double x2) { return pow(x1, 2) - pow(x2, 2) - 1; } double df1_dx1(double x1, double x2) { return 2 * x1; } double df1_dx2(double x1, double x2) { return 2 * x2; } double df2_dx1(double x1, double x2) { return 2 * x1; } double df2_dx2(double x1, double x2) { return -2 * x2; } void newton(double &x1, double &x2) { double eps = 1e-8; int maxIter = 1000; int iter = 0; while (iter < maxIter) { double J[2][2]; J[0][0] = df1_dx1(x1, x2); J[0][1] = df1_dx2(x1, x2); J[1][0] = df2_dx1(x1, x2); J[1][1] = df2_dx2(x1, x2); double F[2]; F[0] = f1(x1, x2); F[1] = f2(x1, x2); double detJ = J[0][0] * J[1][1] - J[0][1] * J[1][0]; double invJ[2][2]; invJ[0][0] = J[1][1] / detJ; invJ[0][1] = -J[0][1] / detJ; invJ[1][0] = -J[1][0] / detJ; invJ[1][1] = J[0][0] / detJ; double dx1 = invJ[0][0] * F[0] + invJ[0][1] * F[1]; double dx2 = invJ[1][0] * F[0] + invJ[1][1] * F[1]; x1 -= dx1; x2 -= dx2; double err = sqrt(dx1 * dx1 + dx2 * dx2); if (err < eps) { break; } iter++; } cout << "x1 = " << x1 << endl; cout << "x2 = " << x2 << endl; } int main() { double x1 = 1.0; double x2 = 1.0; newton(x1, x2); return 0; } ``` 在上述代码中,我们定义了两个线性函数f1和f2,并且分别计算了它们对x1和x2的偏导数。在newton函数中,我们实现牛顿迭代法的迭代公式,并且将计算结果输出。最后,在main函数中,我们初始化了x1和x2,并且调用newton函数求解方程组的解。

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值