We consider PDE constrained nonparametric regression problems in which the parameter $f$ is the unknown coefficient function of a second order elliptic partial differential operator $L_f$, and the unique solution $u_f$ of the boundary value problem \[L_fu=g_1\text{ on } \mathcal O, \quad u=g_2 \text{ on }\partial \mathcal O,\] is observed corrupted by additive Gaussian white noise. Here $\mathcal O$ is a bounded domain in $\mathbb R^d$ with smooth boundary $\partial \mathcal O$, and $g_1, g_2$ are given functions defined on $\mathcal O, \partial \mathcal O$, respectively. Concrete examples include $L_fu=Δu-2fu$ (Schrödinger equation with attenuation potential $f$) and $L_fu=\text{div} (f\nabla u)$ (divergence form equation with conductivity $f$). In both cases, the parameter space \[\mathcal F=\{f\in H^α(\mathcal O)| f > 0\}, ~α>0, \] where $H^α(\mathcal O)$ is the usual order $α$ Sobolev space, induces a set of non-linearly constrained regression functions $\{u_f: f \in \mathcal F\}$. We study Tikhonov-type penalised least squares estimators $\hat f$ for $f$. The penalty functionals are of squared Sobolev-norm type and thus $\hat f$ can also be interpreted as a Bayesian `MAP'-estimator corresponding to some Gaussian process prior. We derive rates of convergence of $\hat f$ and of $u_{\hat f}$, to $f, u_f$, respectively. We prove that the rates obtained are minimax-optimal in prediction loss. Our bounds are derived from a general convergence rate result for non-linear inverse problems whose forward map satisfies a modulus of continuity condition, a result of independent interest that is applicable also to linear inverse problems, illustrated in an example with the Radon transform.