求弹塑性时程分析程序（MATLAB）

hub115

如题
徐赵东《matlab语言在建筑工程中的应用中》一书中油一段弹塑性时程分析程序，本人条了很久也没有调通，不知道有没有人调通过的，发出来大家借鉴一下！

hub115

急切希望大家帮忙看看，这个东西很多人能用到，调通以后对我是个帮助，对其他想用的人也是帮助。

xzning123

楼主你好，我也看了这本书，不过弹塑性那段的程序没有看懂，不知楼主可否注解一下，我倒是能把他的程序调出来，下面给出他的 两个程序
楼主可将我的程序中的elcen.dat换成你的EI centro.dat这样就可以运行了

1. 弹性%%%%%%%%%%%%% MAIN PROGRAM %%%%%%%%%%%%%
   
2. % response analysis of structure--- elastic time history analysis method
   
3. % structure parameter
   
4. tic;
   
5. h=[4000,3300,3300]; % 各层层高, mm
   
6. m=[2.762 2.760 2.300]*1e+3;
   
7. k0=[2.485 1.921 1.522]*1e+5;
   
8. cn=length(m); % 总层数
   
9. % earthquake parameter
   
10. load elcen.dat
    
11. dzhbo=elcen(:,2);   %dzhbo  指地震波
    
12. ct=1.4; % Wilson-θ法的θ值
    
13. dt=0.01;
    
14. xs=200/max(abs(dzhbo)); % 调整地震输入加速度幅值
    
15. ag=dzhbo*0.01*xs; % 地震波, 400 步, 步长为0.02s, 单位m/s2
    
16. ndzh=2000;
    
17. ag1=ag(1:ndzh);
    
18. ag2=ag(2:ndzh+1);
    
19. agtao=ct*(ag2-ag1);
    
20. % initial value
    
21. chsh=zeros(cn,1);
    
22. wyi1=chsh;   %位移
    
23. sdu1=chsh;   %速度
    
24. jsdu1=chsh;  %加速度
    
25. wyimt=chsh;  
    
26. sdumt=chsh;
    
27. jsdumt=chsh;
    
28. unit=ones(cn,1);
    
29. m=diag(m);
    
30. [ik]=matrixju(k0,cn);  %子程序形成刚度矩阵
    
31. [x,d]=eig(ik,m);
    
32. d=sqrt(d);
    
33. w=sort(diag(d));   %将频率从大到小排列
    
34. a=2*w(1)*w(2)*(0.05*w(2)-0.07*w(1))/(w(2)^2-w(1)^2);  %认为是瑞利阻尼，求系数
    
35. b=2*(0.07*w(2)-0.05*w(1))/(w(2)^2-w(1)^2);
    
36. c0=a*m+b*ik;   
    
37. for i=1:ndzh
    
38.     kxin=ik+(3/(ct*dt))*c0+(6/(ct*ct*dt*dt))*m;   %kxin为新的刚度
    
39.     dpxin=-m*unit*agtao(i)+m*(6/(ct*dt)*sdu1+3*jsdu1)+c0*(3*sdu1+ct*dt/2*jsdu1); %新的力增量
    
40.     dxtao=kxin\dpxin;
    
41.     dtjsdu=6*dxtao/(ct*(ct^2*dt^2))-6*sdu1/(ct*ct*dt)-(3/ct)*jsdu1;
    
42.     jsdu=jsdu1+dtjsdu;
    
43.     dtsdu=(dt/2)*(jsdu+jsdu1);
    
44.     sdu=sdu1+dtsdu;
    
45.     dtwyi=dt*sdu1+(1/3)*dt^2*jsdu1+(dt^2/6)*jsdu;
    
46.     wyi=wyi1+dtwyi;
    
47.     jsdu=-m\(m*unit*ag2(i)+c0*sdu+ik*wyi); % 调整加速度
    
48.     wyi1=wyi;
    
49.     sdu1=sdu;
    
50.     jsdu1=jsdu;
    
51.     wyimt=[wyimt wyi*1000];
    
52.     sdumt=[sdumt sdu];
    
53.     jsdumt=[jsdumt jsdu];
    
54. end
    
55. mm=toc;
    
56. mm
    
57. t=0:dt:ndzh*dt;
    
58. figure(2)
    
59. subplot(2,2,1)
    
60. plot(t,wyimt(3,:)/1000,'r-')
    
61. subplot(2,2,2)
    
62. plot(t,jsdumt(3,:),'r-')
    
63. 弹塑性
    
64. %%%%%%%%% MAIN PROGRAM %%%%%%%%%
    
65. %(made in July-October, 2000 by Xu Zhao-dong)
    
66. %(revised in January 6--16, 2002)%
    
67. % this program is three-poly stiffness degeneration model
    
68. % this program is the reinforced concrete structure
    
69. % structure parameters
    
70. m0=[1.690 1.581 1.412]*1e+4;
    
71. k0=[1.512 2.012 2.012]*1e+7;
    
72. h=[4000 7300 10600];
    
73. k3stif=[k0' k0'*0.4 k0'*0.1];
    
74. cn=length(m0);
    
75. m=diag(m0);
    
76. xc=[6.3 4.9 4.2]*1e-3; % 层间开裂位移
    
77. xy=[21.8 18.9 17.2]*1e-3; % 层间屈服位移
    
78. kxi1=0.05; kxi2=0.07; % 阻尼比
    
79. % earthquake parameters
    
80. load elcen.dat
    
81. dzhbo=elcen(:,2);
    
82. ct=1.4;
    
83. dt=0.02;
    
84. ndzh=400;
    
85. tao=ct*dt;
    
86. xs=400/max(abs(dzhbo));
    
87. ag=xs*0.01*dzhbo;
    
88. ag1=ag(1:ndzh);
    
89. ag2=ag(2:ndzh+1);
    
90. agtao=ct*(ag2-ag1);
    
91. % initial value
    
92. chsh=zeros(cn,1);
    
93. pd=chsh;pdd=chsh;
    
94. xcj=chsh;sducj=chsh;jsducj=chsh;fcj=chsh;
    
95. xcjq=chsh;sducjq=chsh;jsducjq=chsh;fcjq=chsh;
    
96. wyi=chsh;sdu=chsh;jsdu=chsh;
    
97. wyimt=chsh;sdumt=chsh;jsdumt=chsh;
    
98. plastic=chsh;plasticmt=chsh;
    
99. unit=ones(cn,1);
    
100. xp=xy;xn=-xy;
     
101. kk=k3stif(:,1);
     
102. [k4stif,fp,fn]=strpara(cn,k3stif,xc,xy); % 求解结构的参数
     
103. for i=1:ndzh
     
104.     pdt=zeros(cn,1);
     
105.     for j=1:cn % 刚度的确定
     
106.     if pd(j)==0
     
107.         kk(j)=k4stif(j,1);
     
108.         plastic(j)=0;
     
109.         fcj(j)=kk(j)*xcj(j)+plastic(j);
     
110.         if xcj(j)>xc(j)
     
111.             pd(j)=2;
     
112.             pdt(j)=1;
     
113.             bl=(xc(j)-xcjq(j,i-1))/(xcj(j)-xcjq(j,i-1));
     
114.             kk(j)=k4stif(j,2);
     
115.             plastic(j)=(k4stif(j,1)-k4stif(j,2))*xc(j);
     
116.             fcj(j)=kk(j)*xcj(j)+plastic(j);
     
117.             elseif xcj(j)<-xc(j)
     
118.                 pd(j)=-2;
     
119.                 pdt(j)=1;
     
120.                 bl=(-xc(j)-xcjq(j,i-1))/(xcj(j)-xcjq(j,i-1));
     
121.                 kk(j)=k4stif(j,2);
     
122.                 plastic(j)=-(k4stif(j,1)-k4stif(j,2))*xc(j);
     
123.                 fcj(j)=kk(j)*xcj(j)+plastic(j);
     
124.         end
     
125.     end
     
126.     if pd(j)==2
     
127.         kk(j)=k4stif(j,2);
     
128.         plastic(j)=(k4stif(j,1)-k4stif(j,2))*xc(j);
     
129.         fcj(j)=kk(j)*xcj(j)+plastic(j);
     
130.         if xcj(j)>xy(j)
     
131.             pd(j)=3;
     
132.             pdt(j)=1;
     
133.             bl=(xy(j)-xcjq(j,i-1))/(xcj(j)-xcjq(j,i-1));
     
134.             kk(j)=k4stif(j,3);
     
135.             plastic(j)=(k4stif(j,4)-k4stif(j,3))*xy(j);
     
136.             fcj(j)=kk(j)*xcj(j)+plastic(j);
     
137.             elseif (xcj(j)<xy(j))&&(sducj(j)<0)
     
138.                 x1a(j)=abs(xcj(j));
     
139.                 pd(j)=1;
     
140.                 kk(j)=((k4stif(j,1)-k4stif(j,2))*xc(j))/x1a(j)+k4stif(j,2);
     
141.                 plastic(j)=0;
     
142.                 fcj(j)=kk(j)*xcj(j)+plastic(j);
     
143.         end
     
144.     end
     
145.     if pd(j)==-2
     
146.         kk(j)=k4stif(j,2);
     
147.         plastic(j)=-(k4stif(j,1)-k4stif(j,2))*xc(j);
     
148.         fcj(j)=kk(j)*xcj(j)+plastic(j);
     
149.         if xcj(j)<-xy(j)
     
150.             pd(j)=-3;
     
151.             pdt(j)=1;
     
152.             bl=(-xy(j)-xcjq(j,i-1))/(xcj(j)-xcjq(j,i-1));
     
153.             kk(j)=k4stif(j,3);
     
154.             plastic(j)=-(k4stif(j,4)-k4stif(j,3))*xy(j);
     
155.             fcj(j)=kk(j)*xcj(j)+plastic(j);
     
156.             elseif (xcj(j)>-xy(j))&&(sducj(j)>0)
     
157.                 x1a(j)=abs(xcj(j));
     
158.                 pd(j)=1;
     
159.                 kk(j)=((k4stif(j,1)-k4stif(j,2))*xc(j))/x1a(j)+k4stif(j,2);
     
160.                 plastic(j)=0;
     
161.                 fcj(j)=kk(j)*xcj(j)+plastic(j);
     
162.         end
     
163.     end
     
164.     if pd(j)==1
     
165.         kk(j)=((k4stif(j,1)-k4stif(j,2))*xc(j))/x1a(j)+k4stif(j,2);
     
166.         plastic(j)=0;
     
167.         fcj(j)=kk(j)*xcj(j)+plastic(j);
     
168.         if xcj(j)>x1a(j)
     
169.             pd(j)=2;
     
170.             pdt(j)=1;
     
171.             bl=(x1a(j)-xcjq(j,i-1))/(xcj(j)-xcjq(j,i-1));
     
172.             kk(j)=k4stif(j,2);
     
173.             plastic(j)=(k4stif(j,1)-k4stif(j,2))*xc(j);
     
174.             fcj(j)=kk(j)*xcj(j)+plastic(j);
     
175.         elseif xcj(j)<-x1a(j)
     
176.             pd(j)=-2;
     
177.             pdt(j)=1;
     
178.             bl=(-x1a(j)-xcjq(j,i-1))/(xcj(j)-xcjq(j,i-1));
     
179.             kk(j)=k4stif(j,2);
     
180.             plastic(j)=-(k4stif(j,1)-k4stif(j,2))*xc(j);
     
181.             fcj(j)=kk(j)*xcj(j)+plastic(j);
     
182.         end
     
183.     end
     
184.     if pd(j)==3
     
185.         kk(j)=k4stif(j,3);
     
186.         plastic(j)=(k4stif(j,4)-k4stif(j,3))*xy(j);
     
187.         fcj(j)=kk(j)*xcj(j)+plastic(j);
     
188.         if sducj(j)<0
     
189.             xp(j)=xcj(j);
     
190.             fp(j)=fcj(j);
     
191.             pd(j)=4;
     
192.             kk(j)=k4stif(j,4);
     
193.             plastic(j)=fp(j)-kk(j)*xp(j);
     
194.             fcj(j)=kk(j)*xcj(j)+plastic(j);
     
195.         end
     
196.     end
     
197.     if pd(j)==-3
     
198.         kk(j)=k4stif(j,3);
     
199.         plastic(j)=-(k4stif(j,4)-k4stif(j,3))*xy(j);
     
200.         fcj(j)=kk(j)*xcj(j)+plastic(j);
     
201.         if sducj(j)>0
     
202.             xn(j)=xcj(j);
     
203.             fn(j)=fcj(j);
     
204.             pd(j)=-4;
     
205.             kk(j)=k4stif(j,4);
     
206.             plastic(j)=fn(j)-kk(j)*xn(j);
     
207.             fcj(j)=kk(j)*xcj(j)+plastic(j);
     
208.         end
     
209.     end
     
210.     if pd(j)==4
     
211.         kk(j)=k4stif(j,4);
     
212.         plastic(j)=fp(j)-kk(j)*xp(j);
     
213.         fcj(j)=kk(j)*xcj(j)+plastic(j);
     
214.         if fcj(j)<0
     
215.             pd(j)=-6;
     
216.             pdt(j)=1;
     
217.             bl=fcjq(j,i-1)/(fcjq(j,i-1)-fcj(j));
     
218.             xa(j)=xcjq(j,i-1)+bl*(xcj(j)-xcjq(j,i-1));
     
219.             fa(j)=0;
     
220.             kk(j)=(fn(j)-fa(j))/(xn(j)-xa(j));
     
221.             plastic(j)=fa(j)-kk(j)*xa(j);
     
222.             fcj(j)=kk(j)*xcj(j)+plastic(j);
     
223.         elseif xcj(j)>xp(j)
     
224.             pd(j)=3;
     
225.             pdt(j)=1;
     
226.             bl=(xp(j)-xcjq(j,i-1))/(xcj(j)-xcjq(j,i-1));
     
227.             kk(j)=k4stif(j,3);
     
228.             plastic(j)=(k4stif(j,4)-k4stif(j,3))*xy(j);
     
229.             fcj(j)=kk(j)*xcj(j)+plastic(j);
     
230.         end
     
231.     end
     
232.     if pd(j)==-4
     
233.         kk(j)=k4stif(j,4);
     
234.         plastic(j)=fn(j)-kk(j)*xn(j);
     
235.         fcj(j)=kk(j)*xcj(j)+plastic(j);
     
236.         if fcj(j)>0
     
237.             pd(j)=6;
     
238.             pdt(j)=1;
     
239.             bl=fcjq(j,i-1)/(fcjq(j,i-1)-fcj(j));
     
240.             xb(j)=xcjq(j,i-1)+bl*(xcj(j)-xcjq(j,i-1));
     
241.             fb(j)=0;
     
242.             kk(j)=(fp(j)-fb(j))/(xp(j)-xb(j));
     
243.             plastic(j)=fb(j)-kk(j)*xb(j);
     
244.             fcj(j)=kk(j)*xcj(j)+plastic(j);
     
245.         elseif xcj(j)<xn(j)
     
246.             pd(j)=-3;
     
247.             pdt(j)=1;
     
248.             bl=(xn(j)-xcjq(j,i-1))/(xcj(j)-xcjq(j,i-1));
     
249.             kk(j)=k4stif(j,3);
     
250.             plastic(j)=-(k4stif(j,4)-k4stif(j,3))*xy(j);
     
251.             fcj(j)=kk(j)*xcj(j)+plastic(j);
     
252.         end
     
253.     end
     
254.     if pd(j)==-6
     
255.         kk(j)=(fn(j)-fa(j))/(xn(j)-xa(j));
     
256.         plastic(j)=fa(j)-kk(j)*xa(j);
     
257.         fcj(j)=kk(j)*xcj(j)+plastic(j);
     
258.         if sducj(j)>0
     
259.             xe(j)=xcj(j);
     
260.             fe(j)=fcj(j);
     
261.             pd(j)=-5;
     
262.             kk(j)=k4stif(j,4);
     
263.             plastic(j)=fe(j)-kk(j)*xe(j);
     
264.             fcj(j)=kk(j)*xcj(j)+plastic(j);
     
265.         elseif xcj(j)<xn(j)
     
266.             pd(j)=-3;
     
267.             pdt(j)=1;
     
268.             bl=(xn(j)-xcjq(j,i-1))/(xcj(j)-xcjq(j,i-1));
     
269.             kk(j)=k4stif(j,3);
     
270.             plastic(j)=-(k4stif(j,4)-k4stif(j,3))*xy(j);
     
271.             fcj(j)=kk(j)*xcj(j)+plastic(j);
     
272.         end
     
273.     end
     
274.     if pd(j)==6
     
275.         kk(j)=(fp(j)-fb(j))/(xp(j)-xb(j));
     
276.         plastic(j)=fb(j)-kk(j)*xb(j);
     
277.         fcj(j)=kk(j)*xcj(j)+plastic(j);
     
278.         if sducj(j)<0
     
279.             xf(j)=xcj(j);
     
280.             ff(j)=fcj(j);
     
281.             pd(j)=5;
     
282.             kk(j)=k4stif(j,4);
     
283.             plastic(j)=ff(j)-kk(j)*xf(j);
     
284.             fcj(j)=kk(j)*xcj(j)+plastic(j);
     
285.         elseif xcj(j)>xp(j)
     
286.             pd(j)=3;
     
287.             pdt(j)=1;
     
288.             bl=(xp(j)-xcjq(j,i-1))/(xcj(j)-xcjq(j,i-1));
     
289.             kk(j)=k4stif(j,3);
     
290.             plastic(j)=(k4stif(j,4)-k4stif(j,3))*xy(j);
     
291.             fcj(j)=kk(j)*xcj(j)+plastic(j);
     
292.         end
     
293.     end  
     
294.     if pd(j)==-5
     
295.         kk(j)=k4stif(j,4);
     
296.         plastic(j)=fe(j)-kk(j)*xe(j);
     
297.         fcj(j)=kk(j)*xcj(j)+plastic(j);
     
298.         if fcj(j)>0
     
299.             pd(j)=6;
     
300.             pdt(j)=1;
     
301.             bl=fcjq(j,i-1)/(fcjq(j,i-1)-fcj(j));
     
302.             xb(j)=xcjq(j,i-1)+bl*(xcj(j)-xcjq(j,i-1));
     
303.             fb(j)=0;
     
304.             kk(j)=(fp(j)-fb(j))/(xp(j)-xb(j));
     
305.             plastic(j)=fb(j)-kk(j)*xb(j);
     
306.             fcj(j)=kk(j)*xcj(j)+plastic(j);
     
307.         elseif (sducj(j)<0)&&(fcj(j)<=0)
     
308.             xa(j)=xcj(j);
     
309.             fa(j)=fcj(j);
     
310.             pd(j)=-6;
     
311.             kk(j)=(fn(j)-fa(j))/(xn(j)-xa(j));
     
312.             plastic(j)=fa(j)-kk(j)*xa(j);
     
313.             fcj(j)=kk(j)*xcj(j)+plastic(j);
     
314.         end
     
315.     end
     
316.     if pd(j)==5
     
317.         kk(j)=k4stif(j,4);
     
318.         plastic(j)=ff(j)-kk(j)*xf(j);
     
319.         fcj(j)=kk(j)*xcj(j)+plastic(j);
     
320.         if fcj(j)<0
     
321.             pd(j)=-6;
     
322.             pdt(j)=1;
     
323.             bl=fcjq(j,i-1)/(fcjq(j,i-1)-fcj(j));
     
324.             xa(j)=xcjq(j,i-1)+bl*(xcj(j)-xcjq(j,i-1));
     
325.             fa(j)=0;
     
326.             kk(j)=(fn(j)-fa(j))/(xn(j)-xa(j));
     
327.             plastic(j)=fa(j)-kk(j)*xa(j);
     
328.             fcj(j)=kk(j)*xcj(j)+plastic(j);
     
329.         elseif (sducj(j)>0)&&(fcj(j)>=0)
     
330.             xb(j)=xcj(j);
     
331.             fb(j)=fcj(j);
     
332.             pd(j)=6;
     
333.             kk(j)=(fp(j)-fb(j))/(xp(j)-xb(j));
     
334.             plastic(j)=fb(j)-kk(j)*xb(j);
     
335.             fcj(j)=kk(j)*xcj(j)+plastic(j);
     
336.         end      
     
337.     end   
     
338.     if pdt(j)==1 % 插值处理
     
339.         dt1=(1-bl)*dt;
     
340.         tao1=ct*dt1;
     
341.         [kju]=matrixju(kk,cn);
     
342.         [c0]=strdamp(m,kju,kxi1,kxi2,cn);
     
343.         dtp1=(1-bl)*dtplastic;
     
344.         for s=1:cn
     
345.             xcjtp(s)=xcjq(s,i-1)+bl*(xcj(s)-xcjq(s,i-1));
     
346.             sducjtp(s)=sducjq(s,i-1)+bl*(sducj(s)-sducjq(s,i-1));
     
347.             jsducjtp(s)=jsducjq(s,i-1)+bl*(jsducj(s)-jsducjq(s,i-1));
     
348.         end     
     
349.         [wyi1]=cjzhuanhuan(xcjtp,cn);
     
350.         [sdu1]=cjzhuanhuan(sducjtp,cn);
     
351.         [jsdu1]=cjzhuanhuan(jsducjtp,cn);
     
352.         kxin1=kju+(3/tao1)*c0+(6/tao1^2)*m;
     
353.         dpxin1=-m*unit*agtao(i-1)*(1-bl)+m*(6/(tao1)*sdu1+3*jsdu1)+c0*(3*sdu1+tao1/2*jsdu1)-ct*dtp1;
     
354.         dxtao1=kxin1\dpxin1;
     
355.         dtjsdu1=6*dxtao1/(ct*(tao1^2))-6*sdu1/(ct*tao1)-(3/ct)*jsdu1;
     
356.         jsdu=jsdu1+dtjsdu1;
     
357.         dtsdu1=(dt1/2)*(jsdu+jsdu1);
     
358.         sdu=sdu1+dtsdu1;
     
359.         dtwyi1=dt1*sdu1+(dt1^2/3)*jsdu1+(dt1^2/6)*jsdu;
     
360.         wyi=wyi1+dtwyi1;
     
361.         [xcj]=zhuanhuan(wyi,cn);
     
362.         [sducj]=zhuanhuan(sdu,cn);
     
363.         [jsducj]=zhuanhuan(jsdu,cn);
     
364.         xcjq(:,i)=xcj;
     
365.         sducjq(:,i)=sducj;
     
366.         jsducjq(:,i)=jsducj;
     
367.         fcjq(:,i)=fcj;
     
368.         wyimt(:,i)=wyi*1000;
     
369.         sdumt(:,i)=sdu;
     
370.         jsdumt(:,i)=jsdu;
     
371.     end   
     
372. end % 1:cn
     
373. pdd=[pdd pd]; % 阶段符号矩阵
     
374. [kju]=matrixju(kk,cn); % 刚度矩阵的聚合
     
375. [c0]=strdamp(m,kju,kxi1,kxi2,cn); % 求解结构的阻尼矩阵
     
376. plastic2=[plastic(2:cn);0];
     
377. plastic3=plastic-plastic2;
     
378. plasticmt=[plasticmt plastic3];
     
379. dtplastic=plastic3-plasticmt(:,i);
     
380. [wyi1]=cjzhuanhuan(xcj,cn);
     
381. [sdu1]=cjzhuanhuan(sducj,cn);
     
382. [jsdu1]=cjzhuanhuan(jsducj,cn);
     
383. kxin=kju+(3/tao)*c0+(6/tao^2)*m;
     
384. dpxin=-m*unit*agtao(i)+m*(6/tao*sdu1+3*jsdu1)+c0*(3*sdu1+tao/2*jsdu1)- ct*dtplastic;
     
385. dxtao=kxin\dpxin;
     
386. dtjsdu=6*dxtao/(ct*(tao^2))-6*sdu1/(ct*tao)-(3/ct)*jsdu1;
     
387. jsdu=jsdu1+dtjsdu; % new acceleration
     
388. dtsdu=(dt/2)*(jsdu+jsdu1);
     
389. sdu=sdu1+dtsdu; % new velocity
     
390. dtwyi=dt*sdu1+(dt^2/3)*jsdu1+(dt^2/6)*jsdu;
     
391. wyi=wyi1+dtwyi; % new displacement
     
392. [xcj]=zhuanhuan(wyi,cn);
     
393. [sducj]=zhuanhuan(sdu,cn);
     
394. jsdu=-m\(m*unit*ag2(i)+c0*sdu+kju*wyi+plastic3); % 加速度修正
     
395. [jsducj]=zhuanhuan(jsdu,cn);
     
396. fcj=kk.*xcj+plastic;
     
397. xcjq=[xcjq xcj];
     
398. sducjq=[sducjq sducj];
     
399. jsducjq=[jsducjq jsducj];
     
400. fcjq=[fcjq fcj];
     
401. wyimt=[wyimt wyi*1000];
     
402. sdumt=[sdumt sdu];
     
403. jsdumt=[jsdumt jsdu];
     
404. end % 1:ndzh
     
405. t=0:0.02:8;
     
406. subplot(2,2,1)
     
407. plot(t,wyimt(3,:),'k');
     
408. subplot(2,2,2)
     
409. plot(t,jsdumt(3,:),'k');

复制代码

xzning123

不好意思，不知道怎么传图片，你运行下就可以看到了

xzning123

可以了，红色为弹性的，下面的为弹塑性的

0911416360

THX FOR YOUR SHARING

fxd3000

想看看，我是菜鸟