#includeiostream#includemath.h/* 设第一条二维曲线方程为 ya*x*xb*xc 则 _ _-1 __ | _ _T _ _ | _ _T _ _ |a| | |x1*x1 x1 1 | |x1*x1 x1 1 | | |x1*x1 x1 1 | | y1 | |b| | |... ... ...| * |... ... ...| | * |... ... ...| * |... | |c| | |xn*xn xn 1 | |xn*xn xn 1 | | |xn*xn xn 1 | | yn | -- | - - - - | - - - - - - 设第二条二维曲线方程为 yd*x*xe*xf 则 _ _-1 __ | _ _T _ _ | _ _T _ _ |d| | |x1*x1 x1 1 | |x1*x1 x1 1 | | |x1*x1 x1 1 | | y1 | |e| | |... ... ...| * |... ... ...| | * |... ... ...| * |... | |f| | |xn*xn xn 1 | |xn*xn xn 1 | | |xn*xn xn 1 | | yn | -- | - - - - | - - - - - - */usingnamespacestd;doublea[][2]{// 第一条曲线周围的点-1,0,0,1,1,0};#defineNsizeof(a)/sizeof(a[0])doubled[][2]{// 第一条曲线周围的点-1,0,0,-1,1,0};#defineNsizeof(a)/sizeof(a[0])#defineMsizeof(d)/sizeof(d[0])voidmatrix_inverse(double(*a)[3],double(*b)[3]);intmain(){// 第一条曲线求系数doublea0[N][3]{0};doublea1[3][N]{0};doublea2[3][3]{0};doublea3[3][3]{0};doublea4[3][N]{0};doubleb0[N][1]{0};doublec0[3][1]{0};for(inti0;iN;i){a0[i][0]pow(a[i][0],2);a0[i][1]a[i][0];a0[i][2]1;a1[0][i]pow(a[i][0],2);// 矩阵转置a1[1][i]a[i][0];a1[2][i]1;b0[i][0]a[i][1];}for(inti0;i3;i)// 矩阵相乘 (k的取值来源于左边矩阵的列数或右边矩阵的行数for(intj0;j3;j)for(intk0;kN;k)a2[i][j](a1[i][k]*a0[k][j]);matrix_inverse(a2,a3);// 矩阵求逆for(inti0;i3;i)for(intj0;jN;j)for(intk0;k3;k)a4[i][j](a3[i][k]*a1[k][j]);for(inti0;i3;i)for(intj0;j1;j)for(intk0;kN;k)c0[i][j](a4[i][k]*b0[k][j]);// for(int i 0; i 3; i)// {// for(int j 0; j 1; j)// {// cout c0[i][j] \t;// }// cout \n;// }// 第二条曲线求系数doubled0[M][3]{0};doubled1[3][M]{0};doubled2[3][3]{0};doubled3[3][3]{0};doubled4[3][M]{0};doublee0[M][1]{0};doublef0[3][1]{0};for(inti0;iM;i){d0[i][0]pow(d[i][0],2);d0[i][1]d[i][0];d0[i][2]1;d1[0][i]pow(d[i][0],2);// 矩阵转置d1[1][i]d[i][0];d1[2][i]1;e0[i][0]d[i][1];}for(inti0;i3;i)// 矩阵相乘 (k的取值来源于左边矩阵的列数或右边矩阵的行数for(intj0;j3;j)for(intk0;kM;k)d2[i][j](d1[i][k]*d0[k][j]);matrix_inverse(d2,d3);// 矩阵求逆for(inti0;i3;i)for(intj0;jM;j)for(intk0;k3;k)d4[i][j](d3[i][k]*d1[k][j]);for(inti0;i3;i)for(intj0;j1;j)for(intk0;kM;k)f0[i][j](d4[i][k]*e0[k][j]);// for(int i 0; i 3; i)// {// for(int j 0; j 1; j)// {// cout f0[i][j] \t;// }// cout \n;// }// 交点所在方程// (c0[0][0] - f0[0][0]) * x*x (c0[0][1] - f0[0][1]) * x (c0[0][2] - f0[0][2]) 0doubleg[2][2]{0};if((pow(c0[0][1]-f0[0][1],2)-4*(c0[0][0]-f0[0][0])*(c0[0][2]-f0[0][2]))0){cout曲线不存在交点\t;return0;}else{g[0][0]-(((c0[0][1]-f0[0][1])pow((pow(c0[0][1]-f0[0][1],2)-4*(c0[0][0]-f0[0][0])*(c0[0][2]-f0[0][2])),0.5))/(2*(c0[0][0]-f0[0][0])));g[0][1](c0[0][0]-f0[0][0])*pow(g[0][0],2)(c0[0][1]-f0[0][1])*g[0][0](c0[0][2]-f0[0][2]);g[1][0]-(((c0[0][1]-f0[0][1])-pow((pow(c0[0][1]-f0[0][1],2)-4*(c0[0][0]-f0[0][0])*(c0[0][2]-f0[0][2])),0.5))/(2*(c0[0][0]-f0[0][0])));g[1][1](c0[0][0]-f0[0][0])*pow(g[1][0],2)(c0[0][1]-f0[0][1])*g[1][0](c0[0][2]-f0[0][2]);}for(inti0;i2;i){for(intj0;j2;j){cout g[i][j]\t;}cout\n;}return0;}voidmatrix_inverse(double(*a)[3],double(*b)[3])// 矩阵求逆{usingnamespacestd;inti,j,k;doublemax,temp;// 定义一个临时矩阵tdoublet[3][3];// 将a矩阵临时存放在矩阵t[n][n]中for(i0;i3;i)for(j0;j3;j)t[i][j]a[i][j];// 初始化B矩阵为单位矩阵for(i0;i3;i)for(j0;j3;j)b[i][j](ij)?(double)1:0;// 进行列主消元找到每一列的主元for(i0;i3;i){maxt[i][i];// 用于记录每一列中的第几个元素为主元ki;// 寻找每一列中的主元元素for(ji1;j3;j){if(fabs(t[j][i])fabs(max)){maxt[j][i];kj;}}//coutthe max number is maxendl;// 如果主元所在的行不是第i行则进行行交换if(k!i){// 进行行交换for(j0;j3;j){tempt[i][j];t[i][j]t[k][j];t[k][j]temp;// 伴随矩阵B也要进行行交换tempb[i][j];b[i][j]b[k][j];b[k][j]temp;}}if(t[i][i]0){cout\nthe matrix does not exist inverse matrix\n;break;}// 获取列主元素tempt[i][i];// 将主元所在的行进行单位化处理//cout\nthe temp is tempendl;for(j0;j3;j){t[i][j]t[i][j]/temp;b[i][j]b[i][j]/temp;}for(j0;j3;j){if(j!i){tempt[j][i];//消去该列的其他元素for(k0;k3;k){t[j][k]t[j][k]-temp*t[i][k];b[j][k]b[j][k]-temp*b[i][k];}}}}}