最简单的Newton-Raphson法求解程序

molen

Newton-Raphson法是有限元非线性分析的基础。
下面的这个FORTRAN程序是最简单的非线性求解演示，一个单自由度的非线性材料桁架用N－R法求解。
很有参考价值。

PROGRAM truss_box_1_2
! -------------------------------------------------------------------------
! Newton-Raphson solver for 1 d.o.f. nonlinear truss example
! Input:
!  d        = horizontal span
!  x_0      = initial height
!  area_0   = initial area
!  E        = Young's modulus
!  n_incr   = number of load increments
!  f_incr   = force increment
!  c_norm   = residual force convergence norm
!  m_iter   = max number of Newton-Raphson iterations
! -------------------------------------------------------------------------
IMPLICIT none
INTEGER,   PARAMETER :: D_P = kind (1.d0) ! machine portable
INTEGER,   PARAMETER :: n_incr = 5, m_iter = 20
REAL(D_P), PARAMETER :: d = 2.5d3, x_0 = 2.5d3, area_0 = 1.d2,     &
                         E = 5.d5, f_incr = 0.3d7, c_norm = 1.d-20, &
                         tolerence = TINY (1.d0)
INTEGER   :: incrm            ! current increment
INTEGER   :: n_iter           ! current N-R iteration number
REAL(D_P) :: force = 0.d0     ! current force
REAL(D_P) :: vertical         ! current vertical truss force
REAL(D_P) :: residual = 0.d0  ! current residual
REAL(D_P) :: r_norm           ! current norm of residual
REAL(D_P) :: length, length_0 ! current and initial length
REAL(D_P) :: stiff            ! current stiffness
REAL(D_P) :: stress           ! current stress
REAL(D_P) :: x, area          ! current height and area
REAL(D_P) :: volume           ! current volume
REAL(D_P) :: u                ! current displacement increment
!  initialize geometry and stiffness
  x = x_0 ; area = area_0
  length_0 = sqrt ( x **2 + d **2 )
  volume   = area_0 * length_0
  stiff    = (area / length_0) * E * (x / length_0 ) **2
!  start load increments
  DO incrm = 1, n_incr
    force    = force + f_incr
    residual = residual - f_incr
!  Newton-Raphson iteration
    r_norm = c_norm * 2 ; n_iter = 0  ! initialize
    DO WHILE ( (r_norm > c_norm) .and. &
               (n_iter < m_iter)       )
      n_iter = n_iter + 1
!     Find geometry increment
      u      = -residual / stiff      ! displacement increment
      x      = x + u                  ! update height
      length = sqrt ( x **2 + d **2 ) ! update length
      area   = volume / length        ! update area
      
!     Find stresses and residual forces
      stress   = E * log (length / length_0)  ! axial stress
      vertical = stress * area * x / length   ! vertical force
      residual = vertical - force
      r_norm   = abs (residual / force)
      print "('increment = ', i3, ' iteration = ', i3,     &
      &       ' residual = ', g10.4, ' x = ', g10.4,       &
      &       ' force = ', g10.4)", incrm, n_iter, r_norm, &
      &       x, force
!     Find stiffness and check stability
      stiff = (area/length) * (E - 2 * stress) * (x/length) **2 &
            + (stress * area / length)
      IF ( abs ( stiff) < tolerence )  THEN
        print *, 'Near zero stiffness'
        STOP     'Near zero stiffness'
      END IF ! unstable
    END DO ! Newton-Raphson iteration
  END DO ! incrm
END PROGRAM truss_box_1_2
!
! Running with parameters E = 5.d5, f_incr = 1.5d7, c_norm = 1.d-20:
! increment = 1 iteration = 1 residual = 0.2184     x = 4621. force = 0.1500E+08
! increment = 1 iteration = 2 residual = 0.4239E-01 x = 5540. force = 0.1500E+08
! increment = 1 iteration = 3 residual = 0.2973E-02 x = 5822. force = 0.1500E+08
! increment = 1 iteration = 4 residual = 0.1806E-04 x = 5844. force = 0.1500E+08
! increment = 1 iteration = 5 residual = 0.6780E-09 x = 5845. force = 0.1500E+08
! increment = 1 iteration = 6 residual = 0.1242E-15 x = 5845. force = 0.1500E+08
! increment = 1 iteration = 7 residual =  0.000     x = 5845. force = 0.1500E+08

[ 本帖最后由 molen 于 2007-12-21 07:28 编辑 ]

bellgreen

什么语言哈？楼主还不如给各伪语言格式的更加的通用

jian4dan

楼主你能用汉语给一点一点解释一下吗？十分感谢

苦雨孤鸿

楼主一开始不就说了嘛 ，是Fortran的

kzhu

3# jian4dan
楼主都把源程序给出了，还叫人家解释，有点过分了，有点Fortran基础的应该都能看懂的

bbssbb

quadratic convergence, NICE!

jsl024110

对自己编程进行非线性计算，彻底掌握Newton-Raphson方法大有裨益，顶版主提供的程序。要是能将弧长法以简单实例展示出来就更完美了！

guayigou

fortran 语言呢，无需解释了吧，玩有限元的，这个必须的呀

syky86

FORTRAN语言不大懂，解释一下就好了

zdqwh

确实很有有参考价值

tjz1981

没有看出来是Fortran语言啊

daniel1985

谢谢楼主分享，让我省事多了

zengyu113

谢谢楼主分享

redplum1

是用Fortran编写的，多谢多谢

fireflywph

看不懂，没有学过Fortran

dabing437

刚度表达式 看不懂啥意思