61 lines
1.6 KiB
C
61 lines
1.6 KiB
C
|
|
extern int mm,n,mpt;
|
|
extern float h,c2,anorm,x[];
|
|
|
|
void difeq(k,k1,k2,jsf,is1,isf,indexv,ne,s,y)
|
|
float **s,**y;
|
|
int indexv[],is1,isf,jsf,k,k1,k2,ne;
|
|
{
|
|
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];
|
|
}
|
|
}
|