difeq.c

适合大型数值计算代码 现在网络上已经找不到了 购买需要20$

extern int mm,n,mpt;
extern float h,c2,anorm,x[];

void difeq(int k, int k1, int k2, int jsf, int is1, int isf, int indexv[],
	int ne, float **s, float **y)
{
	float temp,temp1,temp2;

	if (k == k1) {
		if (n+mm & 1) {
			s[3][3+indexv[1]]=1.0;
			s[3][3+indexv[2]]=0.0;
			s[3][3+indexv[3]]=0.0;
			s[3][jsf]=y[1][1];
		} else {
			s[3][3+indexv[1]]=0.0;
			s[3][3+indexv[2]]=1.0;
			s[3][3+indexv[3]]=0.0;
			s[3][jsf]=y[2][1];
		}
	} else if (k > k2) {
		s[1][3+indexv[1]] = -(y[3][mpt]-c2)/(2.0*(mm+1.0));
		s[1][3+indexv[2]]=1.0;
		s[1][3+indexv[3]] = -y[1][mpt]/(2.0*(mm+1.0));
		s[1][jsf]=y[2][mpt]-(y[3][mpt]-c2)*y[1][mpt]/(2.0*(mm+1.0));
		s[2][3+indexv[1]]=1.0;
		s[2][3+indexv[2]]=0.0;
		s[2][3+indexv[3]]=0.0;
		s[2][jsf]=y[1][mpt]-anorm;
	} else {
		s[1][indexv[1]] = -1.0;
		s[1][indexv[2]] = -0.5*h;
		s[1][indexv[3]]=0.0;
		s[1][3+indexv[1]]=1.0;
		s[1][3+indexv[2]] = -0.5*h;
		s[1][3+indexv[3]]=0.0;
		temp1=x[k]+x[k-1];
		temp=h/(1.0-temp1*temp1*0.25);
		temp2=0.5*(y[3][k]+y[3][k-1])-c2*0.25*temp1*temp1;
		s[2][indexv[1]]=temp*temp2*0.5;
		s[2][indexv[2]] = -1.0-0.5*temp*(mm+1.0)*temp1;
		s[2][indexv[3]]=0.25*temp*(y[1][k]+y[1][k-1]);
		s[2][3+indexv[1]]=s[2][indexv[1]];
		s[2][3+indexv[2]]=2.0+s[2][indexv[2]];
		s[2][3+indexv[3]]=s[2][indexv[3]];
		s[3][indexv[1]]=0.0;
		s[3][indexv[2]]=0.0;
		s[3][indexv[3]] = -1.0;
		s[3][3+indexv[1]]=0.0;
		s[3][3+indexv[2]]=0.0;
		s[3][3+indexv[3]]=1.0;
		s[1][jsf]=y[1][k]-y[1][k-1]-0.5*h*(y[2][k]+y[2][k-1]);
		s[2][jsf]=y[2][k]-y[2][k-1]-temp*((x[k]+x[k-1])
			*0.5*(mm+1.0)*(y[2][k]+y[2][k-1])-temp2
			*0.5*(y[1][k]+y[1][k-1]));
		s[3][jsf]=y[3][k]-y[3][k-1];
	}
}