请教Newmark方法C++程序？谢谢！

tianya7071 · 发表于 2007-7-27 10:25

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

您需要登录才可以下载或查看，没有账号？我要加入

x

谁有Newmark方法C++程序，做课题急用，请大虾指点一下，邮箱tianyarrrbbb@163.com

风花雪月 · 发表于 2007-7-29 15:32

http://sokocalo.engr.ucdavis.edu ... penSees_source_.tgz

[ 本帖最后由风花雪月于 2007-7-29 15:35 编辑 ]

hao1982 · 发表于 2007-8-3 17:40

#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*[n];
for(i=0;i<n;i++)
M[i]=new double[n];
double **C=new double*[n];
for(i=0;i<n;i++)
C[i]=new double[n];
double **K=new double*[n];
for(i=0;i<n;i++)
K[i]=new double[n];
double **X=new double*[nt+1];
for(i=0;i<nt+1;i++)
X[i]=new double[n];
double **V=new double*[nt+1];
for(i=0;i<nt+1;i++)
V[i]=new double[n];
double **A=new double*[nt+1];
for(i=0;i<nt+1;i++)
A[i]=new double[n];
double *P=new double[n];
X[0][0]=0;X[0][1]=0;
V[0][0]=0;V[0][1]=0;
A[0][0]=0;A[0][1]=2;
M[0][0]=1;M[0][1]=0;M[1][0]=0;M[1][1]=1;
C[0][0]=0;C[0][1]=0;C[1][0]=0;C[1][1]=0;
K[0][0]=0;K[0][1]=-1;K[1][0]=1;K[1][1]=0;
P[0]=0;P[1]=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[i][j]<<setw(12)<<V[i][j]<<setw(12)<<A[i][j]<<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[size];
double *temp=new double[size];
double *tran=new double[size];
double **KK=new double*[size];
for(i=0;i<size;i++)
KK[i]=new double[size];
double **L=new double* [size];
for(i=0;i<size;i++)
L[i]=new double [size];
double **U=new double* [size];
for(i=0;i<size;i++)
U[i]=new double [size];
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[i][j]=K[i][j]+a0*M[i][j]+a1*C[i][j];
}
ArrayLU(KK,L,U,size);//LU Reduce
for(ti=1;ti<=nt;ti++)
{
for(i=0;i<size;i++)
temp[i]=a0*X[ti-1][i]+a2*V[ti-1][i]+a3*A[ti-1][i];
ArrayMVector(M,temp,tran,size);
for(i=0;i<size;i++)
PT[i]=P[i]+tran[i];
for(i=0;i<size;i++)
temp[i]=a1*X[ti-1][i]+a4*V[ti-1][i]+a5*A[ti-1][i];
ArrayMVector(C,temp,tran,size);
for(i=0;i<size;i++)
PT[i]=PT[i]+tran[i];
LUSolve(L,U,PT,temp,size);
for(i=0;i<size;i++)
{
X[ti][i]=temp[i];
A[ti][i]=a0*(X[ti][i]-X[ti-1][i])-a2*V[ti-1][i]-a3*A[ti-1][i];
V[ti][i]=V[ti-1][i]+a6*A[ti-1][i]+a7*A[ti][i];
}
}
for(i=0;i<size;i++)
delete [] KK[i];
delete []KK;
for(i=0;i<size;i++)
delete [] L[i];
delete []L;
for(i=0;i<size;i++)
delete [] U[i];
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[size];
Y[0]=P[0]/L[0][0];
for(i=1;i<size;i++)
{
sum=0;
for(j=0;j<i;j++)
sum+=L[i][j]*Y[j];
Y[i]=(P[i]-sum)/L[i][i];
}
X[size-1]=Y[size-1]/U[size-1][size-1];
for(i=size-2;i>=0;i--)
{
sum=0;
for(j=i+1;j<size;j++)
sum+=U[i][j]*X[j];
X[i]=(Y[i]-sum)/U[i][i];
}
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* [size];
for(i=0;i<size;i++)
B[i]=new double [size];
for(i=0;i<size;i++)
{
for(j=0;j<size;j++)
{
B[i][j]=0;L[i][j]=0;U[i][j]=0;
}
B[i][i]=1;
}
for(i=0;i<size-1;i++)
{
if(A[i][i]==0)
{
cout<<"the Array is singular"<<endl;
}
for(j=i+1;j<size;j++)
{
scale=A[j][i]/A[i][i];
for(k=i;k<size;k++)
A[j][k]=A[j][k]-A[i][k]*scale;
for(k=0;k<size;k++)
B[j][k]=B[j][k]-B[i][k]*scale;
}
}
Inverse(B,L,size);
for(i=0;i<size;i++)
{
for(j=i;j<size;j++)
U[i][j]=A[i][j];
}
for(i=0;i<size;i++)
delete [] B[i];
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[i][i]=1;
for(i=0;i<size-1;i++)
{
if(A[i][i]==0)
{
cout<<"the Array is singular"<<endl;
}
for(j=i+1;j<size;j++)
{
scale=A[j][i]/A[i][i];
for(k=i;k<size;k++)
A[j][k]=A[j][k]-A[i][k]*scale;
for(k=0;k<size;k++)
B[j][k]=B[j][k]-B[i][k]*scale;
}
}
for(i=size-1;i>0;i--)
{
for(k=0;k<size;k++)
B[i][k]=B[i][k]/A[i][i];
A[i][i]=1;
for(j=i-1;j>=0;j--)
{
scale=A[j][i]/A[i][i];
A[j][i]=0;
for(k=0;k<size;k++)
B[j][k]=B[j][k]-B[i][k]*scale;
}
}
for(k=0;k<size;k++)
B[0][k]=B[0][k]/A[0][0];
}
//***************************************************************************************************
//***************************************************************************************************
//***************************************************************************************************
//****************** Array Multiply Vector ******************
//***************************************************************************************************
void ArrayMVector(double **A,double *X,double *B,int size)
{
for(int i=0;i<size;i++)
{
B[i]=0;
for(int j=0;j<size;j++)
{
B[i]+=A[i][j]*X[j];
}
}
}
//***************************************************************************************************
//***************************************************************************************************

复制代码

[ 本帖最后由风花雪月于 2007-9-3 09:51 编辑 ]

mort · 发表于 2011-11-27 11:09

有具体详细简单的这方面的程序吗？希望高人指点

账号		自动登录	找回密码
密码			我要加入

[经典算法] 请教Newmark方法C++程序？谢谢！

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

本帖被以下淘专辑推荐:

评分