Abstract: A·U+B(U)=f, wherein A is a linear differential operator, B is a nonlinear differential operator, and f is an inhomogeneous term (source term) in a nonlinear partial differential equation to be satisfied by a physical quantity U, is solved by successive approximation. In calculation, (f−A·Um−B(Um)) is given as a nonlinear residual rr of an approximate solution Um, wherein m is the number of repeating times, and the approximate solution Um is repeatedly corrected so as to reduce a norm of a nonlinear residual rm+1 employed in a subsequent step.