This paper addresses the effects of non-linear, non-Darcy flows in reservoir on the value of the productivity index. The productivity index (PI) of the well draining a reservoir with no flux on the exterior boundaries is defined as the ratio between the production rate and the pressure drawdown. Experience shows that during the dynamical process of hydrocarbon recovery this ratio stabilizes to some constant value, which, in general, is a non-linear function of both the pressure drawdown and the production rate. The linear (Darcy) case is well understood, and excellent approximate formulae are available to compute the PI for various well-reservoir geometries. These formulae are generally obtained through a semi-analytical solution of the transient problem. To handle the more realistic non-linear situation, the current practice is to solve first the nonlinear problem many times for different values of production rate, and then to add some ad-hoc corrective parameter(s) in the linear formulae in order to reproduce the non-linear nature of the flow. Our approach, based on recent progress in the modeling of transient Forchheimer flows, uses partly non-numerical techniques to evaluate the productivity index. It provides, for a wide class of reservoir geometries, an accurate enough approximate algebraic relation between the PI for the non-linear Darcy-Forchheimer flows and the PI for the Darcy flows. We show that the solution of the original problem can be obtained by applying direct variational methods, which are computationally less expensive than the grid based techniques. In addition, by using the existing and newly developed approaches in the theory of symmetrization, alternative "non-analytical" algorithms are presented to assist optimal well design, without repeatedly solving numerically the reservoir simulation problem. The approach will be demonstrated on practical well optimization problems.