Методы решения жестких и нежестких краевых задач

for(nn=1;nn<=nn_last;nn++){ //A CYCLE BY HARMONICS STARTING FROM THE 1st HARMONIC (EXCEPT ZERO ONE)

double x=0.0;//A coordinate from the left boundary - it is needed in case of heterogeneous system of the ODE for computing the particular vector FF

double expo_from_minus_step[8][8]={0};//The matrix for placement of the exponent in it at the step of (0-x1) type

double expo_from_plus_step[8][8]={0};// The matrix for placement of the exponent in it in the step of (x1-0) type

double mat_row_for_minus_expo[8][8]={0};//the auxiliary matrix for particular vector computing when moving at step of (0-x1) type

double mat_row_for_plus_expo[8][8]={0};// the auxiliary matrix for particular vector computing when moving at step of (x1-0) type

double U[4][8]={0};//The matrix of the boundary conditions of the left boundary of the dimensionality 4x8

double u_[4]={0};//Dimensionality 4 vector of the external influence for the boundary conditions of the left boundary

double V[4][8]={0};//The boundary conditions matrix of the right boundary of the dimensionality 4x8

double v_[4]={0};// The dimensionality 4 vector of the external influence for the boundary conditions for the right boundary

double Y[100+1][8]={0};//The array of the vector-solutions of the corresponding linear algebraic equations system (in each point of the interval between the boundaries): MATRIXS*Y=VECTORS

double A[8][8]={0};//Matrix of the coefficients of the system of ODE

double FF[8]={0};//Vector of the particular solution of the heterogeneous ODE at the integration interval sector

double Y_many[8*11]={0};// a composite vector consisting of the vectors Y(xi) in 11 points from point 0 (left boundary Y(0) to the point 10 (right boundary Y(x10))

double MATRIX_many[8*11][8*11]={0};//The matrix of the system of the linear algebraic equations

double B_many[8*11]={0};// a vector of the right parts of the SLAE: MATRIX_many*Y_many=B_many

double Y_vspom[8]={0};//an auxiliary vector

double Y_rezult[8]={0};//an auxiliary vector

double nn2,nn3,nn4,nn5,nn6,nn7,nn8;//Number of the nn harmonic raised to corresponding powers

nn2=nn*nn; nn3=nn2*nn; nn4=nn2*nn2; nn5=nn4*nn; nn6=nn4*nn2; nn7=nn6*nn; nn8=nn4*nn4;

//Filling-in of non-zero elements of the A matrix of the coefficients of ODE system

A[0][1]=1.0;

A[1][0]=(1-nju)/2*nn2; A[1][3]=-(1+nju)/2*nn; A[1][5]=-nju;

A[2][3]=1.0;

A[3][1]=(1+nju)/(1-nju)*nn; A[3][2]=2*nn2/(1-nju); A[3][4]=2*nn/(1-nju);

A[4][5]=1.0;

A[5][6]=1.0;

A[6][7]=1.0;

A[7][1]=-nju/c2; A[7][2]=-nn/c2; A[7][4]=-(nn4+1/c2); A[7][6]=2*nn2;

//Here firstly it is necessary to make filling-in by non-zero values of the matrix and boundary conditions vector U*Y[0]=u_ (on the left) and V*Y[100]=v_ (on the right):

U[0][1]=1.0; U[0][2]=nn*nju; U[0][4]=nju; u_[0]=0.0;//Force T1 at the left boundary is equal to zero

U[1][0]=-(1-nju)/2*nn; U[1][3]=(1-nju)/2; U[1][5]=(1-nju)*nn*c2; u_[1]=0.0;//Force S* at the left boundary is equal to zero

U[2][4]=-nju*nn2; U[2][6]=1.0; u_[2]=0;//Momentum M1 at the left boundary is equal to zero

U[3][5]=(2-nju)*nn2; U[3][7]=-1.0;

u_[3]=-sin(nn*gamma)/(nn*gamma);//Force Q1* at the left boundary is distributed at the angle -gamma +gamma

V[0][0]=1.0; v_[0]=0.0;// The right boundary displacement u is equal to zero

V[1][2]=1.0; v_[1]=0.0;//The right boundary displacement v is equal to zero

V[2][4]=1.0; v_[2]=0.0;//The right boundary displacement w is equal to zero

V[3][5]=1.0; v_[3]=0.0;//The right boundary rotation angle is equal to zero

//Here initial filling-in of U*Y[0]=u_ и V*Y[100]=v_ terminates

exponent(A,(-step*10),expo_from_minus_step);//A negative step (step value is less than zero due to direction of matrix exponent computation)

//x=0.0;//the initial value of coordinate - for partial vector computation

//mat_row_for_partial_vector(A, step, mat_row_for_minus_expo);

//Filling-in of the SLAE coefficients matrix MATRIX_many

for(int i=0;i<4;i++){

for(int j=0;j<8;j++){

MATRIX_many[i][j]=U[i][j];

MATRIX_many[8*11-4+i][8*11-8+j]=V[i][j];

}

B_many[i]=u_[i];

B_many[8*11-4+i]=v_[i];

}

for(int kk=0;kk<(11-1);kk++){//(11-1) unit matrix and EXPO matrix should be written into MATRIX_many

for(int i=0;i<8;i++){

MATRIX_many[i+4+kk*8][i+kk*8]=1.0;//filling-in by unit matrixes

for(int j=0;j<8;j++){

MATRIX_many[i+4+kk*8][j+8+kk*8]=-expo_from_minus_step[i][j];//filling-in by matrix exponents

}

}

}

//Solution of the system of linear algebraic equations

GAUSS(MATRIX_many,B_many,Y_many);

//Computation of the state vectors in 101 points - left point 0 and right point 100 exponent(A,step,expo_from_plus_step);

for(int i=0;i<11;i++){//passing points filling-in in all 10 intervals (we will obtain points from 0 to 100) between 11 nodes

for(int j=0;j<8;j++){

Y[0+i*10][j]=Y_many[j+i*8];//in 11 nodes the vectors are taken from SLAE solutions - from Y many

}

}

for(int i=0;i<10;i++){//passing points filling-in in 10 intervals

for(int j=0;j<8;j++){

Y_vspom[j]=Y[0+i*10][j];//the initial vector for ith segment, zero point, starting point of the ith segment

}

mat_on_vect(expo_from_plus_step, Y_vspom, Y_rezult);

for(int j=0;j<8;j++){

Y[0+i*10+1][j]=Y_rezult[j];//the first point of the interval filling-in

Y_vspom[j]=Y_rezult[j];//for the next step

}

mat_on_vect(expo_from_plus_step, Y_vspom, Y_rezult);

for(int j=0;j<8;j++){

Y[0+i*10+2][j]=Y_rezult[j];//filling-in of the 2nd point of the interval

Y_vspom[j]=Y_rezult[j];//for the next step

}

mat_on_vect(expo_from_plus_step, Y_vspom, Y_rezult);

for(int j=0;j<8;j++){

Y[0+i*10+3][j]=Y_rezult[j];//filling-in of the 3rd point of the interval

Y_vspom[j]=Y_rezult[j];//for the next step

}

mat_on_vect(expo_from_plus_step, Y_vspom, Y_rezult);

for(int j=0;j<8;j++){

Y[0+i*10+4][j]=Y_rezult[j];//filing-in of the 4th point of the interval

Y_vspom[j]=Y_rezult[j];//for the next step

}

mat_on_vect(expo_from_plus_step, Y_vspom, Y_rezult);

for(int j=0;j<8;j++){

Y[0+i*10+5][j]=Y_rezult[j];// filing-in of the 5th point of the interval

Y_vspom[j]=Y_rezult[j];// for the next step

}

mat_on_vect(expo_from_plus_step, Y_vspom, Y_rezult);

for(int j=0;j<8;j++){

Y[0+i*10+6][j]=Y_rezult[j];// filing-in of the 6th point of the interval

Y_vspom[j]=Y_rezult[j];// for the next step

}

mat_on_vect(expo_from_plus_step, Y_vspom, Y_rezult);

for(int j=0;j<8;j++){

Y[0+i*10+7][j]=Y_rezult[j];// filing-in of the 7th point of the interval

Y_vspom[j]=Y_rezult[j];// for the next step

}

mat_on_vect(expo_from_plus_step, Y_vspom, Y_rezult);

for(int j=0;j<8;j++){

Y[0+i*10+8][j]=Y_rezult[j];// filing-in of the 8th point of the interval

Y_vspom[j]=Y_rezult[j];// for the next step

}

mat_on_vect(expo_from_plus_step, Y_vspom, Y_rezult);

for(int j=0;j<8;j++){

Y[0+i*10+9][j]=Y_rezult[j];// filing-in of the 9th point of the interval

Y_vspom[j]=Y_rezult[j];// for the next step

}

}

//Computation of the momentum in all points between the boundaries

for(int ii=0;ii<=100;ii++){

Moment[ii]+=Y[ii][4]*(-nju*nn2)+Y[ii][6]*1.0;//Momentum M1 in the point [ii]

//U[2][4]=-nju*nn2; U[2][6]=1.0; u_[2]=0;//Momentum M1

}

}//CYCLE BY HARMONICS TERMINATES HERE

for(int ii=0;ii<=100;ii++){

fprintf(fp,"%f\n",Moment[ii]);

}

fclose(fp);

printf( "PRESS any key to continue...\n" );

_getch();

return 0;

}

My web-site related to stiff boundary value problems solution methods:

Размещено на Allbest.ru