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

tianya7071

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

风花雪月

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

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

hao1982

1. #include<iostream.h>
   
2. #include<iomanip.h>
   
3. 
   
4. void ArrayMVector(double**,double*,double*,int);
   
5. void ArrayLU(double**,double**,double**,int);
   
6. void Inverse(double**,double**,int);
   
7. void LUSolve(double**,double**,double*,double*,int);
   
8. void Newmark(double**,double**,double**,double*,double**,double**,double**,int,int,double);
   
9. 
   
10. void main()
    
11. {
    
12.         int i,j,n,nt;
    
13.     n=2;nt=100;
    
14.         double dt=0.08;
    
15.         double **M=new double*[n];
    
16.         for(i=0;i<n;i++)
    
17.                 M[i]=new double[n];
    
18.     double **C=new double*[n];
    
19.         for(i=0;i<n;i++)
    
20.                 C[i]=new double[n];
    
21.     double **K=new double*[n];
    
22.         for(i=0;i<n;i++)
    
23.                 K[i]=new double[n];
    
24.     double **X=new double*[nt+1];
    
25.         for(i=0;i<nt+1;i++)
    
26.                 X[i]=new double[n];
    
27.     double **V=new double*[nt+1];
    
28.         for(i=0;i<nt+1;i++)
    
29.                 V[i]=new double[n];
    
30.     double **A=new double*[nt+1];
    
31.         for(i=0;i<nt+1;i++)
    
32.                 A[i]=new double[n];
    
33.         double *P=new double[n];
    
34.         X[0][0]=0;X[0][1]=0;
    
35.         V[0][0]=0;V[0][1]=0;
    
36.         A[0][0]=0;A[0][1]=2;
    
37.         M[0][0]=1;M[0][1]=0;M[1][0]=0;M[1][1]=1;
    
38.     C[0][0]=0;C[0][1]=0;C[1][0]=0;C[1][1]=0;
    
39.     K[0][0]=0;K[0][1]=-1;K[1][0]=1;K[1][1]=0;
    
40.         P[0]=0;P[1]=2;
    
41.     Newmark(M,C,K,P,X,V,A,n,nt,dt);
    
42.         for(i=0;i<=nt;i++)
    
43.         {
    
44.         cout<<i*dt<<" s:"<<endl;
    
45.                 for(j=0;j<n;j++)
    
46.                 cout<<setw(12)<<X[i][j]<<setw(12)<<V[i][j]<<setw(12)<<A[i][j]<<endl;
    
47.         }
    
48. }
    
49. 
    
50. 
    
51. //***************************************************************************************************
    
52. //*************                                Newmark                                  *************
    
53. //***************************************************************************************************
    
54. void Newmark(double**M,double**C,double**K,double *P,double**X,double**V,double**A,int size,int nt,double dt)
    
55. {
    
56.         int ti,i,j;
    
57.     double a=0.25,d=0.5;
    
58.         double a0,a1,a2,a3,a4,a5,a6,a7;
    
59.         double *PT=new double[size];
    
60.     double *temp=new double[size];
    
61.     double *tran=new double[size];
    
62.     double **KK=new double*[size];
    
63.         for(i=0;i<size;i++)
    
64.                 KK[i]=new double[size];
    
65.     double **L=new double* [size];
    
66.         for(i=0;i<size;i++)
    
67.                 L[i]=new double [size];
    
68.     double **U=new double* [size];
    
69.         for(i=0;i<size;i++)
    
70.                 U[i]=new double [size];
    
71. 
    
72.     a0=1/(a*dt*dt);
    
73.         a1=d/(a*dt);
    
74.         a2=1/(a*dt);
    
75.         a3=1/(2*a)-1;
    
76.         a4=d/a-1;
    
77.         a5=dt*(d/a-2)/2;
    
78.         a6=dt*(1-d);
    
79.         a7=d*dt;
    
80.         for(i=0;i<size;i++)
    
81.         {
    
82.                 for(j=0;j<size;j++)
    
83.                         KK[i][j]=K[i][j]+a0*M[i][j]+a1*C[i][j];
    
84.         }
    
85.     ArrayLU(KK,L,U,size);//LU Reduce
    
86.         for(ti=1;ti<=nt;ti++)
    
87.         {
    
88.                 for(i=0;i<size;i++)
    
89.                         temp[i]=a0*X[ti-1][i]+a2*V[ti-1][i]+a3*A[ti-1][i];
    
90.         ArrayMVector(M,temp,tran,size);
    
91.                 for(i=0;i<size;i++)
    
92.                         PT[i]=P[i]+tran[i];
    
93.         for(i=0;i<size;i++)
    
94.                         temp[i]=a1*X[ti-1][i]+a4*V[ti-1][i]+a5*A[ti-1][i];
    
95.         ArrayMVector(C,temp,tran,size);
    
96.                 for(i=0;i<size;i++)
    
97.                         PT[i]=PT[i]+tran[i];
    
98.                 LUSolve(L,U,PT,temp,size);
    
99.                 for(i=0;i<size;i++)
    
100.                 {
     
101.                         X[ti][i]=temp[i];
     
102.                         A[ti][i]=a0*(X[ti][i]-X[ti-1][i])-a2*V[ti-1][i]-a3*A[ti-1][i];
     
103.                         V[ti][i]=V[ti-1][i]+a6*A[ti-1][i]+a7*A[ti][i];
     
104.                 }
     
105.         }
     
106.         for(i=0;i<size;i++)
     
107.                 delete [] KK[i];
     
108.         delete []KK;
     
109.     for(i=0;i<size;i++)
     
110.                 delete [] L[i];
     
111.         delete []L;
     
112.     for(i=0;i<size;i++)
     
113.                 delete [] U[i];
     
114.         delete []U;
     
115.         delete []PT;
     
116.         delete []temp;
     
117.         delete []tran;
     
118. }
     
119. 
     
120. 
     
121. 
     
122. //***************************************************************************************************
     
123. //***************************************************************************************************
     
124. 
     
125. //***************************************************************************************************
     
126. //*************                            LU Reduce of Array                           *************
     
127. //***************************************************************************************************
     
128. void LUSolve(double**L,double**U,double*P,double*X,int size)
     
129. {
     
130.         int i,j;
     
131.         double sum;
     
132.         double* Y=new double[size];
     
133.     Y[0]=P[0]/L[0][0];
     
134.     for(i=1;i<size;i++)
     
135.         {
     
136.         sum=0;
     
137.                 for(j=0;j<i;j++)
     
138.                         sum+=L[i][j]*Y[j];
     
139.                 Y[i]=(P[i]-sum)/L[i][i];
     
140.         }
     
141.     X[size-1]=Y[size-1]/U[size-1][size-1];
     
142.         for(i=size-2;i>=0;i--)
     
143.         {
     
144.                 sum=0;
     
145.                 for(j=i+1;j<size;j++)
     
146.                         sum+=U[i][j]*X[j];
     
147.                 X[i]=(Y[i]-sum)/U[i][i];
     
148.         }
     
149.         delete []Y;
     
150. }
     
151. 
     
152. //***************************************************************************************************
     
153. //***************************************************************************************************
     
154. 
     
155. 
     
156. 
     
157. 
     
158. 
     
159. 
     
160. //***************************************************************************************************
     
161. //*************                            LU Reduce of Array                           *************
     
162. //***************************************************************************************************
     
163. void ArrayLU(double**A,double**L,double**U,int size)
     
164. {
     
165.     int i,j,k;
     
166.         double scale;
     
167.     double **B=new double* [size];
     
168.         for(i=0;i<size;i++)
     
169.                 B[i]=new double [size];
     
170.     for(i=0;i<size;i++)
     
171.         {
     
172.                 for(j=0;j<size;j++)
     
173.                 {
     
174.                         B[i][j]=0;L[i][j]=0;U[i][j]=0;
     
175.                 }
     
176.         B[i][i]=1;
     
177.         }
     
178.         for(i=0;i<size-1;i++)
     
179.         {
     
180.                 if(A[i][i]==0)
     
181.                 {
     
182.                         cout<<"the Array is singular"<<endl;
     
183.                 }
     
184.                 for(j=i+1;j<size;j++)
     
185.                 {
     
186.                         scale=A[j][i]/A[i][i];
     
187.                         for(k=i;k<size;k++)
     
188.                                 A[j][k]=A[j][k]-A[i][k]*scale;
     
189.                         for(k=0;k<size;k++)
     
190.                                 B[j][k]=B[j][k]-B[i][k]*scale;
     
191.                 }
     
192.         }
     
193.         Inverse(B,L,size);
     
194.         for(i=0;i<size;i++)
     
195.         {
     
196.                 for(j=i;j<size;j++)
     
197.                         U[i][j]=A[i][j];
     
198.         }
     
199.         for(i=0;i<size;i++)
     
200.                 delete [] B[i];
     
201.         delete []B;
     
202. }
     
203. //***************************************************************************************************
     
204. //***************************************************************************************************
     
205. 
     
206. 
     
207. 
     
208. 
     
209. //***************************************************************************************************
     
210. //******************                         Array Inverse                         ******************
     
211. //***************************************************************************************************
     
212. void Inverse(double**A,double**B,int size)
     
213. {
     
214.         int i,j,k;
     
215.         double scale;
     
216.         for(i=0;i<size;i++)
     
217.                 B[i][i]=1;
     
218.         for(i=0;i<size-1;i++)
     
219.         {
     
220.                 if(A[i][i]==0)
     
221.                 {
     
222.                         cout<<"the Array is singular"<<endl;
     
223.                 }
     
224.                 for(j=i+1;j<size;j++)
     
225.                 {
     
226.                         scale=A[j][i]/A[i][i];
     
227.                         for(k=i;k<size;k++)
     
228.                                 A[j][k]=A[j][k]-A[i][k]*scale;
     
229.                         for(k=0;k<size;k++)
     
230.                                 B[j][k]=B[j][k]-B[i][k]*scale;
     
231.                 }
     
232.         }
     
233.         for(i=size-1;i>0;i--)
     
234.         {
     
235.                 for(k=0;k<size;k++)
     
236.                         B[i][k]=B[i][k]/A[i][i];
     
237.                 A[i][i]=1;
     
238.                 for(j=i-1;j>=0;j--)
     
239.                 {
     
240.                         scale=A[j][i]/A[i][i];
     
241.                         A[j][i]=0;
     
242.             for(k=0;k<size;k++)
     
243.                                 B[j][k]=B[j][k]-B[i][k]*scale;
     
244.                 }
     
245.         }
     
246.         for(k=0;k<size;k++)
     
247.                 B[0][k]=B[0][k]/A[0][0];
     
248. }
     
249. //***************************************************************************************************
     
250. //***************************************************************************************************
     
251. 
     
252. 
     
253. 
     
254. 
     
255. //***************************************************************************************************
     
256. //******************                      Array Multiply Vector                    ******************
     
257. //***************************************************************************************************
     
258. void ArrayMVector(double **A,double *X,double *B,int size)
     
259. {
     
260.         for(int i=0;i<size;i++)
     
261.         {
     
262.                 B[i]=0;
     
263.                 for(int j=0;j<size;j++)
     
264.                 {
     
265.                         B[i]+=A[i][j]*X[j];
     
266.                 }
     
267.         }
     
268. }
     
269. //***************************************************************************************************
     
270. //***************************************************************************************************

复制代码

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

mort

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