#include <math.h>
#include<iostream>
using namespace std;
extern double *function_rung(int n1,double t1,double *u,double v)
{
double **matrixA =NULL;
double **matrixM=NULL;
double **matrixM_1=NULL;
double **matrixK=NULL;
double **matrixC =NULL;
double **matrixK1=NULL;
double **matrixC1 =NULL;
double *vector_B=NULL;
double *vector_f=NULL;
double *vector_f1=NULL;
double *w1=NULL;
double *u1=NULL;
double *v_k=NULL;
double c1=11593.1014;
double k1=280000.0;
double k2=156000.0;
double m1=12000.0;
double m2=500.0;
double m=6067.0;
double mb=182010.0;
double pi=3.1415926;
double l=30.0;
double temp=0.0;
double EI=2658069000;
double x=0.0;
double x1=0.0;
matrixA=new double *[2*n1+4];
for(int i=0;i<2*n1+4;i++)
{
matrixA[i]=new double [2*n1+4];
}
matrixM=new double *[n1+2];
for(int i=0;i<n1+2;i++)
{
matrixM[i]=new double [n1+2];
}
matrixM_1=new double *[n1+2];
for(int i=0;i<n1+2;i++)
{
matrixM_1[i]=new double [2*n1+4];
}
matrixK=new double *[n1+2];
for(int i=0;i<n1+2;i++)
{
matrixK[i]=new double [n1+2];
}
matrixC=new double *[n1+2];
for(int i=0;i<n1+2;i++)
{
matrixC[i]=new double [n1+2];
}
matrixK1=new double *[n1+2];
for(int i=0;i<n1+2;i++)
{
matrixK1[i]=new double [n1+2];
}
matrixC1=new double *[n1+2];
for(int i=0;i<n1+2;i++)
{
matrixC1[i]=new double [n1+2];
}
w1=new double [n1];
vector_f= new double [n1+2];
vector_f1= new double [n1+2];
u1=new double [2*n1+4];
vector_B= new double [2*n1+4];
v_k=new double [2*n1+4];
// give the value to the matrix M
for(int i=0;i<n1;i++)
{
for(int j=0;j<n1;j++)
{
if(i!=j) matrixM[i][j]=0.0;
matrixM[i][i]=1;
}
}
for(int i=n1;i<n1+2;i++)
{
for(int j=0;j<n1+2;j++)
{
matrixM[i][j]=0.0;
}
}
for(int i=0;i<n1;i++)
{
matrixM[i][n1]=2*(m1/mb)*sin((i+1)*pi*v*t1/30);
}
for(int i=0;i<n1;i++)
{
matrixM[i][n1+1]=sin((i+1)*pi*v*t1/30)*2*(m2/mb);
}
matrixM[n1][n1]=m1;
matrixM[n1+1][n1+1]=m2;
// give the value to the matrix K
for(int i=0;i<n1;i++)
{
w1[i]=(i+1)*(i+1)*pi*pi*sqrt(EI/m)/(l*l);
}
for(int i=0;i<n1;i++)
{
for(int j=0;j<n1+2;j++)
{
if(i!=j) matrixK[i][j]=0;
matrixK[i][i]=w1[i]*w1[i];
}
}
for(int i=0;i<n1;i++)
{
matrixK[n1][i]=0;
}
matrixK[n1][n1]=k1;
matrixK[n1][n1+1]=-k1;
for(int i=0;i<n1;i++)
{
matrixK[n1+1][i]=-1*k2*sin((i+1)*pi*v*t1/l) ;
}
matrixK[n1+1][n1]=-k1;
matrixK[n1+1][n1+1]=k1+k2;
//give the value to matrix C
for(int i=0;i<n1+2;i++)
{
for(int j=0;j<n1+2; j++)
{
matrixC[i][j]=0;
}
}
matrixC[n1][n1]=c1;
matrixC[n1][n1+1]=-c1;
matrixC[n1+1][n1+1]=c1;
matrixC[n1+1][n1]=-c1;
//give the value to vector f
for(int i=0;i<n1;i++)
{
vector_f[i]=(-1)*(2/m)*(m1+m2)*10*sin((i+1)*pi*v*t1/l);
}
vector_f[n1]=0;
vector_f[n1+1]=0;
// solve the ob-matrix
{
double temp= 0;
for(int i=0;i<n1+2;i++)
{
for(int j=0;j<2*n1+4;j++)
{
if (j<n1+2) matrixM_1[i][j]=matrixM[i][j];
else if (j==2+i+n1) matrixM_1[i][j]=1;
else matrixM_1[i][j]=0;
}
}
for(int i=0;i<n1+2;i++)
{
for (int j=0;j<2*n1+4;j++)
{
for(int k=i+1;k<n1+2;k++)
{
matrixM_1[i][j]=matrixM_1[i][j]+matrixM_1[k][j];
}
}
temp=matrixM_1[i][i];
for(int j=0;j<2*n1+4;j++)
{
matrixM_1[i][j]/=temp;
}
for(int k=0;k<n1+2;k++)
{
if(k!=i)
{
temp=matrixM_1[k][i];
for(int l=0;l<2*n1+4;l++)
{
matrixM_1[k][l]=matrixM_1[k][l]-matrixM_1[i][l]*temp;
}
}
}
}
for(int i=0;i<n1+2;i++)
{
for(int j=0;j<2*n1+4;j++)
{
matrixM_1[i][j]=matrixM_1[i][j]/matrixM_1[i][i];
}
}
for(int i=0;i<n1+2;i++)
{
for(int j=0;j<n1+2;j++)
{
matrixM[i][j]=matrixM_1[i][j+n1+2];
}
}
}
// compute the multiplication between matrix
for(int i=0;i<n1+2;i++)
{
for(int j=0;j<n1+2;j++)
{
for(int k=0;k<n1+2;k++)
{
x+=matrixM[i][k]*matrixK[k][j];
}
matrixK1[i][j]=x;
x=0;
}
}
for(int i=0;i<n1+2;i++)
{
for(int j=0;j<n1+2;j++)
{
for(int k=0;k<n1+2;k++)
{
x+=matrixM[i][k]*matrixC[k][j];
}
matrixC1[i][j]=x;
x=0;
}
}
//compute the multiplication between matrix and vector
for(int i=0;i<n1+2;i++)
{
for(int j=0;j<n1+2;j++)
{
x1+=matrixM[i][j]*vector_f[j];
}
vector_f1[i]=x1;
x1=0;
}
// give the value to matrix A
for(int i=0;i<n1+2;i++)
{
for(int j=0;j<n1+2;j++)
{
matrixA[i][j]=0;
}
for(int j=n1+2;j<2*n1+4;j++)
{
matrixA[i][i+n1+2]=1;
if(i!=j-n1-2)matrixA[i][j]=0;
}
}
for(int i=n1+2;i<2*n1+4;i++)
{
for(int j=0;j<n1+2;j++)
{
matrixA[i][j]=-1*matrixK1[i-n1-2][j];
}
}
for(int i=n1+2;i<2*n1+4;i++)
{
for(int j=n1+2;j<2*n1+4;j++)
{
matrixA[i][j]=-1*matrixC1[i-n1-2][j-n1-2];
}
}
// give the value to vector B
for(int i=0;i<n1+2;i++)
{
vector_B[i]=0;
}
for(int i=n1+2;i<2*n1+4;i++)
{
vector_B[i]=vector_f1[i-n1-2];
}
//compute the value
for(int i=0;i<2*n1+4;i++)
{
for(int j=0;j<2*n1+4;j++)
{
x+=matrixA[i][j]*u[j];
}
u1[i]=x;
x=0;
u1[i]+=vector_B[i];
v_k[i]=u1[i];
}
return v_k;
// freedom the space of points
for(int i=0;i<2*n1+4;i++)
{
delete matrixA[i];
}
delete matrixA;
for(int i=0;i<n1+2;i++)
{
delete matrixM[i];
}
delete matrixM;
for(int i=0;i<n1+2;i++)
{
delete matrixM_1[i];
}
delete matrixM_1;
for(int i=0;i<n1+2;i++)
{
delete matrixK[i];
}
delete matrixK;
for(int i=0;i<n1+2;i++)
{
delete matrixC[i];
}
delete matrixC;
for(int i=0;i<n1+2;i++)
{
delete matrixK1[i];
}
delete matrixK1;
for(int i=0;i<n1+2;i++)
{
delete matrixC1[i];
}
delete matrixC1;
delete vector_B;
delete vector_f;
delete vector_f1;
delete u1;
delete v_k;
delete w1;
}
