请教Newmark方法C++程序？谢谢！
谁有Newmark方法C++程序，做课题急用，请大虾指点一下，邮箱tianyarrrbbb@163.com http://sokocalo.engr.ucdavis.edu/~jeremic/OpenSees/SOURCE/_OpenSees_source_.tgz[ 本帖最后由 风花雪月 于 2007-7-29 15:35 编辑 ] #include<iostream.h>
#include<iomanip.h>
void ArrayMVector(double**,double*,double*,int);
void ArrayLU(double**,double**,double**,int);
void Inverse(double**,double**,int);
void LUSolve(double**,double**,double*,double*,int);
void Newmark(double**,double**,double**,double*,double**,double**,double**,int,int,double);
void main()
{
int i,j,n,nt;
n=2;nt=100;
double dt=0.08;
double **M=new double*;
for(i=0;i<n;i++)
M=new double;
double **C=new double*;
for(i=0;i<n;i++)
C=new double;
double **K=new double*;
for(i=0;i<n;i++)
K=new double;
double **X=new double*;
for(i=0;i<nt+1;i++)
X=new double;
double **V=new double*;
for(i=0;i<nt+1;i++)
V=new double;
double **A=new double*;
for(i=0;i<nt+1;i++)
A=new double;
double *P=new double;
X=0;X=0;
V=0;V=0;
A=0;A=2;
M=1;M=0;M=0;M=1;
C=0;C=0;C=0;C=0;
K=0;K=-1;K=1;K=0;
P=0;P=2;
Newmark(M,C,K,P,X,V,A,n,nt,dt);
for(i=0;i<=nt;i++)
{
cout<<i*dt<<" s:"<<endl;
for(j=0;j<n;j++)
cout<<setw(12)<<X<<setw(12)<<V<<setw(12)<<A<<endl;
}
}
//***************************************************************************************************
//************* Newmark *************
//***************************************************************************************************
void Newmark(double**M,double**C,double**K,double *P,double**X,double**V,double**A,int size,int nt,double dt)
{
int ti,i,j;
double a=0.25,d=0.5;
double a0,a1,a2,a3,a4,a5,a6,a7;
double *PT=new double;
double *temp=new double;
double *tran=new double;
double **KK=new double*;
for(i=0;i<size;i++)
KK=new double;
double **L=new double* ;
for(i=0;i<size;i++)
L=new double ;
double **U=new double* ;
for(i=0;i<size;i++)
U=new double ;
a0=1/(a*dt*dt);
a1=d/(a*dt);
a2=1/(a*dt);
a3=1/(2*a)-1;
a4=d/a-1;
a5=dt*(d/a-2)/2;
a6=dt*(1-d);
a7=d*dt;
for(i=0;i<size;i++)
{
for(j=0;j<size;j++)
KK=K+a0*M+a1*C;
}
ArrayLU(KK,L,U,size);//LU Reduce
for(ti=1;ti<=nt;ti++)
{
for(i=0;i<size;i++)
temp=a0*X+a2*V+a3*A;
ArrayMVector(M,temp,tran,size);
for(i=0;i<size;i++)
PT=P+tran;
for(i=0;i<size;i++)
temp=a1*X+a4*V+a5*A;
ArrayMVector(C,temp,tran,size);
for(i=0;i<size;i++)
PT=PT+tran;
LUSolve(L,U,PT,temp,size);
for(i=0;i<size;i++)
{
X=temp;
A=a0*(X-X)-a2*V-a3*A;
V=V+a6*A+a7*A;
}
}
for(i=0;i<size;i++)
delete [] KK;
delete []KK;
for(i=0;i<size;i++)
delete [] L;
delete []L;
for(i=0;i<size;i++)
delete [] U;
delete []U;
delete []PT;
delete []temp;
delete []tran;
}
//***************************************************************************************************
//***************************************************************************************************
//***************************************************************************************************
//************* LU Reduce of Array *************
//***************************************************************************************************
void LUSolve(double**L,double**U,double*P,double*X,int size)
{
int i,j;
double sum;
double* Y=new double;
Y=P/L;
for(i=1;i<size;i++)
{
sum=0;
for(j=0;j<i;j++)
sum+=L*Y;
Y=(P-sum)/L;
}
X=Y/U;
for(i=size-2;i>=0;i--)
{
sum=0;
for(j=i+1;j<size;j++)
sum+=U*X;
X=(Y-sum)/U;
}
delete []Y;
}
//***************************************************************************************************
//***************************************************************************************************
//***************************************************************************************************
//************* LU Reduce of Array *************
//***************************************************************************************************
void ArrayLU(double**A,double**L,double**U,int size)
{
int i,j,k;
double scale;
double **B=new double* ;
for(i=0;i<size;i++)
B=new double ;
for(i=0;i<size;i++)
{
for(j=0;j<size;j++)
{
B=0;L=0;U=0;
}
B=1;
}
for(i=0;i<size-1;i++)
{
if(A==0)
{
cout<<"the Array is singular"<<endl;
}
for(j=i+1;j<size;j++)
{
scale=A/A;
for(k=i;k<size;k++)
A=A-A*scale;
for(k=0;k<size;k++)
B=B-B*scale;
}
}
Inverse(B,L,size);
for(i=0;i<size;i++)
{
for(j=i;j<size;j++)
U=A;
}
for(i=0;i<size;i++)
delete [] B;
delete []B;
}
//***************************************************************************************************
//***************************************************************************************************
//***************************************************************************************************
//****************** Array Inverse ******************
//***************************************************************************************************
void Inverse(double**A,double**B,int size)
{
int i,j,k;
double scale;
for(i=0;i<size;i++)
B=1;
for(i=0;i<size-1;i++)
{
if(A==0)
{
cout<<"the Array is singular"<<endl;
}
for(j=i+1;j<size;j++)
{
scale=A/A;
for(k=i;k<size;k++)
A=A-A*scale;
for(k=0;k<size;k++)
B=B-B*scale;
}
}
for(i=size-1;i>0;i--)
{
for(k=0;k<size;k++)
B=B/A;
A=1;
for(j=i-1;j>=0;j--)
{
scale=A/A;
A=0;
for(k=0;k<size;k++)
B=B-B*scale;
}
}
for(k=0;k<size;k++)
B=B/A;
}
//***************************************************************************************************
//***************************************************************************************************
//***************************************************************************************************
//****************** Array Multiply Vector ******************
//***************************************************************************************************
void ArrayMVector(double **A,double *X,double *B,int size)
{
for(int i=0;i<size;i++)
{
B=0;
for(int j=0;j<size;j++)
{
B+=A*X;
}
}
}
//***************************************************************************************************
//***************************************************************************************************
[ 本帖最后由 风花雪月 于 2007-9-3 09:51 编辑 ] 有具体详细简单的这方面的程序吗？希望高人指点
页:
[1]