ニューマーク-β法のプログラム(サブルーチン)


      subroutine NMB(m,c,k,u,ud,n,na,Ts,Te,dt,beta,iout)

      real m,k
      dimension m(na,na),c(na,na),k(na,na),u(na),ud(na)
	dimension f(50),u2d(50),a(50,50),p(50),q(50),r(50)
	nb=50
 
c---- calculation of u2d at time Ts

	do i=1,n
	do j=1,n
		a(i,j)=m(i,j)
	enddo
	enddo

	call fwr(a,n,nb)

	call excite(f,n,nb,Ts)

	do i=1,n
		r(i)=f(i)
		do j=1,n
			r(i)=r(i)-c(i,j)*ud(j)-k(i,j)*u(j)
		enddo
	enddo

	call bws(a,r,n,nb)

	do i=1,n
		u2d(i)=r(i)
	enddo

c	write( 6,*) Ts,u(iout),ud(iout),u2d(iout)
	write(10,*) Ts,u(iout),ud(iout),u2d(iout)

c---- set constants ---
	c1=dt/2.
	c2=beta*dt**2
	c3=(0.5-beta)*dt**2
	c4=dt**2/2.
	c5=beta*dt**2

	do i=1,n
	do j=1,n
		a(i,j)=m(i,j)+c1*c(i,j)+c2*k(i,j)
	enddo
	enddo
	call fwr(a,n,nb)

c---- clock ---
      itmax=ifix((Te-Ts)/dt+0.1)+1

	do it=2,itmax
		t=Ts+float(it-1)*dt

c---- calculation of u2d at time t
		do i=1,n
			p(i)=ud(i)+c1*u2d(i)
			q(i)=u(i)+ud(i)*dt+c3*u2d(i)
		enddo

		call excite(f,n,nb,t)

		do i=1,n
			r(i)=f(i)
			do j=1,n
				r(i)=r(i)-c(i,j)*p(j)-k(i,j)*q(j)
			enddo
		enddo

		call bws(a,r,n,nb)
		
c---- calculation of u,ud at time t

		do i=1,n
			u2d(i)=r(i)
			u(i)=u(i)+ud(i)*dt+c4*u2d(i)+c5*(r(i)-u2d(i))
			ud(i)=ud(i)+c1*(r(i)+u2d(i))
		enddo
		
c		write( 6,*) t,u(iout),ud(iout),u2d(iout)
		write(10,*) t,u(iout),ud(iout),u2d(iout)
		
	enddo
	
      return
	end



	subroutine fwr(a,n,na)

	dimension a(na,1)

	do k=1,n-1
		do j=k+1,n
			a(k,j)=a(k,j)/a(k,k)
		enddo

		do i=k+1,n
		do j=k+1,n
			a(i,j)=a(i,j)-a(i,k)*a(k,j)
		enddo
		enddo
	enddo

	return
	end



	subroutine bws(a,x,n,na)

	dimension a(na,1),x(na)

	do k=1,n-1
		x(k)=x(k)/a(k,k)
		do i=k+1,n
			x(i)=x(i)-a(i,k)*x(k)
		enddo
	enddo
	
     	x(n)=x(n)/a(n,n)

	do l=1,n-1
		k=n-l
		do j=k+1,n
			x(k)=x(k)-a(k,j)*x(j)
		enddo
	enddo

	return
	end