关于newmark法问题的请教！
subroutine newmark(count,WMM,WDM,WKM,WEF,U,DiffU,DiffDiffU)use constant
use crank
use imsl
implicit none
integer i,j !下标
real(8) WMM(61,61), WDM(61,61), WKM(61,61),WEF(61) ！输入
real(8) U(0:400,61),DiffU(0:400,61),DiffDiffU(0:400,61) ！输出位移、速度、加速度
real(8) tempU(1:11,61),tempDiffU(1:11,61),tempDiffDiffU(1:11,61) ！中间变量：位移、速度、加速度
real(8) temp1(61),temp2(61),TempWMM(61,61),TempWDM(61,61) ！中间变量
real(8) K(61,61),P(61) ！中间变量
integer count ！代表所求时刻0.1、0.2、0.3**********
real(8)::a=5,b=0.25,dt=0.01
real(8) A0,A1,A2,A3,A4,A5,A6,A7,crankrad,temp ！参数
A0=1/(b*dt**2)
A1=a/(b*dt)
A2=1/(b*dt)
A3=1/(2*b)-1
A4=a/b-1
A5=dt*a/(2*b)-dt
A6=(1-a)*dt
A7=a*dt
crankrad=crankspeed*2*pi/60
temp=crankrad*180/(9*pi)
WDM=WDM*temp ！进行自变量的转换，与方法本身无关
WMM=WMM*temp**2 ！进行自变量的转换，与方法本身无关
tempU(1,:)=U(count-1,:) ！设定初始值
tempDiffU(1,:)=DiffU(count-1,:)
tempDiffDiffU(1,:)=DiffDiffU(count-1,:)
TempWMM=A0*WMM
TempWDM=A1*WDM
k=TempWMM+TempWDM+WKM
do i=2,11
temp1=A0*tempU(i-1,:)+A2*tempDiffU(i-1,:)+A3*tempDiffDiffU(i-1,:)
temp2=A1*tempU(i-1,:)+A4*tempDiffU(i-1,:)+A5*tempDiffDiffU(i-1,:)
temp1=WMM.x.temp1
temp2=WDM.x.temp2
p=WEF+temp1+temp2
tempU(i,:)=k.ix.p
tempDiffDiffU(i,:)=A0*(tempU(i,:)-tempU(i-1,:))-A2*tempDiffU(i-1,:)-A3*tempDiffDiffU(i-1,:)
tempDiffU(i,:)=tempU(i-1,:)+A6*tempDiffDiffU(i-1,:)+A7*tempDiffDiffU(i,:)
enddo
U(count,:)=tempU(11,:)
DiffU(count,:)=tempDiffU(11,:)
DiffDiffU(count,:)=tempDiffDiffU(11,:)
do i=1,11
write(8,"(61f60.8)")tempU(i,:)
enddo
return
end subroutine newmark 我这个是在fortran中编的。
问题是：不收敛
说明：
（1）我无论怎么调我的时间步长，都没有用；反而是我时间步长越小，发散的越快
（2）我无论怎么改变参数a、b都不起效果
（3）我初始条件在第一次调用子函数时是自己任意给的，但也试了零和非零的没有效果。后面再调用时就是前一次计算的结果
最后恳请前辈们能帮我分析一下，到底是哪里出了问题？
万分感激！ 看看建立的力学模型是否有问题。
[ 本帖最后由 无水1324 于 2007-11-27 13:54 编辑 ] 刚度矩阵、质量矩阵、阻尼矩阵我都检查了，没检查出毛病。
关键是程序中的temp1，temp2变化很大，刚开始只有0.000001这样子，后来增加到10000000数量级，后面还有更大的，不可思议级了。
我这个程序编的对吗？
这个程序的含义是这样的：
我要对某个机构进行机构弹性动力学分析，现在呢就是想求出机构运转一周各个杆件的弹性变形过程。
我把一个周期分成40个等分，在某一个等分里假设构件的刚度矩阵、质量矩阵、阻尼矩阵都是不变的；再对每一等分应用newmark方法求解，初始值即为前一等分的最后的值，这样是否正确。
还有就是我直接对整个周期应用newmark方法，请问这两种应用有什么区别吗，我最后算出来的结果是不一样的，虽然都是发散的。
[ 本帖最后由 spring_zhao 于 2007-11-27 14:01 编辑 ] 建议先用一个简单的例子调试，待方法本身通过以后，再带入复杂的模型运算。这样便于检查。 哦，这个建议不错，我这就去试一试 Program newmark
implicit none
integer i
real(8) j !下标
real(8) WMM, WDM, WKM,WEF
real(8) U(0:30),DiffU(0:30),DiffDiffU(0:30)
real(8) K,P
real(8)::a=5,b=0.25,dt=0.1
real(8) A0,A1,A2,A3,A4,A5,A6,A7
A0=1/(b*dt**2)
A1=a/(b*dt)
A2=1/(b*dt)
A3=1/(2*b)-1
A4=a/b-1
A5=dt*a/(2*b)-dt
A6=(1-a)*dt
A7=a*dt
U(0)=0
DiffU(0)=0
DiffDiffU(0)=0
WMM=0.02
WDM=0.028
do i=1,30
j=i*0.10000000000000000
if(i>=0.and.i<=10) then
WKM=1
else
WKM=0.1
endif
WEF=5*(1-j/2.0)**2
k=A0*WMM+A1*WDM+WKM
p=WEF+WMM*(A0*U(i-1)+A2*DiffU(i-1)+A3*DiffDiffU(i-1))+&
WDM*(A1*U(i-1)+A4*DiffU(i-1)+A5*DiffDiffU(i-1))
U(i)=p/k
DiffDiffU(i)=A0*(U(i)-U(i-1))-A2*DiffU(i-1)-A3*DiffDiffU(i-1)
DiffU(i)=U(i-1)+A6*DiffDiffU(i-1)+A7*DiffDiffU(i)
write(*,"('time=',3f12.5)")j,WEF,U(i)
enddo
end Program newmark
这是我的简单例子，但是算出来的结果和正确的结果相差很大，应该是不正确的，请问增量的和我给出的这种newmaik法有什么区别吗，我倒是觉得没什么区别 DiffDiffU(i)=A0*(U(i)-U(i-1))-A2*DiffU(i-1)-A3*DiffDiffU(i-1)
DiffU(i)=U(i-1)+A6*DiffDiffU(i-1)+A7*DiffDiffU(i)
我发现数值算法，这两个应变量的不同写法，会带来不同的结果 可能是我的算法出了问题，我发现同样的非线性系统，用全量的newmark法计算是发散的，而且怎么调b都没有用，而用增量的newmark法计算时就可以收敛。
到底是什么机理导致他们这样，还有待于研究。
[ 本帖最后由 spring_zhao 于 2007-11-28 11:19 编辑 ] 终于找到问题了，全量remark法，也是收敛的。
但是有一个问题，全量和增量有什么区别呢，请问前辈们，他们的区别表现在哪呢？ 原来还是算法的问题
谢谢wanyeqing2003 看看这个就是我得出的结果，表示杆件上某点的振动，不知道是不是收敛了，它的含义说明什么？
请高手过来看看阿，都没有人和我讨论问题？:'( 回复 10 # spring_zhao 的帖子
你好，我与你遇到了同样的问题， 能说说你是怎么解决的问题吗？是算法的问题吗还是编程出的问题？谢谢回复！！！！ 楼主你好，我在用newmark法解非线性微分方程组，里面恢复力项有一个非线性回滞环，用newmark法时估计需要用到迭代，求楼主指教！
页:
[1]
2