fuction.cpp

提供矩阵以及向量的各种运算

#pragma once
#include<iostream>
#include"Matrix.h"
using std::cout;
using std::cin;
using std::endl;

double insert(Matrix& A,double x,double y)
{
     if(A.getRow()==4&&A.getLine()==4)
	 {
		 double dx=A.MatrA[0][2]-A.MatrA[0][1];
		 double dy=A.MatrA[2][0]-A.MatrA[1][0];
		 int ki;
		 int kj;
		 double res=0;
		 for(int i=1;i<4;i++)
		 {
			 ki = (i==2)?(-2):1;
			 for(int j=1;j<4;j++)
			 {
				 kj = (j==2)?(-2):1;
				 res+=ki*kj*(y-A.MatrA[A.plus(i,1)][0])*(y-A.MatrA[A.plus(i,2)][0])*(x-A.MatrA[0][A.plus(j,1)])
					 *(x-A.MatrA[0][A.plus(j,2)])*A.MatrA[i][j]/(4*dx*dx*dy*dy);

			 }
		 }
		 return res;
	 }
	 else 
		 cout<< "error!"<<endl;
}



Matrix& descompose(Matrix& A,double x,double y)
{
	if(A.getRow()==7&&A.getLine()==7)
	{
		Matrix* res = new Matrix(4,4);
		int ki,kj;
		int ii=1,jj=1;
		do
		{
			jj++;
			if(jj==A.getLine()-2)
				break;
		}while((x-A.MatrA[0][jj])>fabs(x-A.MatrA[0][jj+1]));
		do
		{
			ii++;
			if(ii==A.getRow()-2)
				break;
		}while((y-A.MatrA[ii][0])>fabs(y-A.MatrA[ii+1][0]));


		res->MatrA[0][0]=0;
		for(int i=1;i<4;i++)
		{
			ki = i==2?0:(i==1?(-1):1);
			res->MatrA[0][i]=A.MatrA[0][jj+ki];
			res->MatrA[i][0]=A.MatrA[ii+ki][0];
			for(int j=1;j<4;j++)
			{

				kj = j==2?0:(j==1?(-1):1);
				res->MatrA[i][j]=A.MatrA[ii+ki][jj+kj];
			}
		}
		//cout<<ii<<endl;
		//cout<<jj<<endl;
		return *res;
	}
	else 
		cout<<"error"<<endl;
}
double compute(Matrix& A,double x,double y)
{
	double res;
	
	res=insert(descompose(A,x,y),x,y);
	return res;
}



vector& nonlinear(double x,double y)
{
	vector X1(4);
	vector X2(4);
	Matrix dF(4,4);
	vector F(4);
	for(int i=0;i<4;i++)
	{
		for(int j=0;j<4;j++)
		{
			dF.MatrA[i][j]=1;
		}
	}
	dF.MatrA[2][0]=0.5;
	dF.MatrA[3][1]=0.5;
	double t=-1;
	int k=1;
	X1.setVector(0,0.328341);
	X1.setVector(1,1.342025);
	X1.setVector(2,-0.502958);
	X1.setVector(3,1.457645);
	do
	{
		X2=X2.copy(X1);
		F.setVector(0,0.5*cos(X1.getVector(0))+X1.getVector(1)+X1.getVector(2)+X1.getVector(3)-2.67-x);
		F.setVector(1,X1.getVector(0)+0.5*sin(X1.getVector(1))+X1.getVector(2)+X1.getVector(3)-1.07-y);
		F.setVector(2,0.5*X1.getVector(0)+X1.getVector(1)+cos(X1.getVector(2))+X1.getVector(3)-3.74-x);
		F.setVector(3,X1.getVector(0)+0.5*X1.getVector(1)+X1.getVector(2)+sin(X1.getVector(3))-0.79-y);
		dF.MatrA[0][0]=-0.5*sin(X1.getVector(0));
		dF.MatrA[1][1]=0.5*cos(X1.getVector(1));
		dF.MatrA[2][2]=-sin(X1.getVector(2));
		dF.MatrA[3][3]=cos(X1.getVector(3));
		X1=X1+(!(dF)*F)/t;
		k++;
		
		if(k==30)
		{
			cout<<"error"<<endl;
			break;
		}
	}while((X1+X2/t).max()/X1.max()>1e-12);
	//cout<<k<<endl;
	vector* res=new vector;
	res->copy(X1);
	return* res;
}


Matrix& approach(Matrix& A,vector& X,vector& Y,int k)
{
	if(X.getN()==A.row&&Y.getN()==A.line)
	{
		Matrix B(A.row,k+1);
		Matrix G(A.line,k+1);
		for(int i=0;i<A.row;i++)
		{
			for(int j=0;j<k+1;j++)
			{
				B.MatrA[i][j]=pow(X.getVector(i),j);
			}
		}
		for(int i=0;i<A.line;i++)
		{
			for(int j=0;j<k+1;j++)
			{
				G.MatrA[i][j]=pow(Y.getVector(i),j);
			}
		}
		//B.output();
		//G.output();
		Matrix C;
		C=(!((~B)*B))*(~B)*A*G*(!((~G)*G));
		Matrix* res=new Matrix;
		res->copy(C);
		return *res;
	}
	else
		cout<<"error"<<endl;
}

double appcompute(Matrix& C,double x,double y)
{
	double res=0;
	for(int i=0;i<C.row;i++)
	{
		
		for(int j=0;j<C.line;j++)
		{
			res+=C.MatrA[i][j]*pow(x,i)*pow(y,j);
		}
	}
	return res;
}