Finite difference for nonlinear PDE

I'm trying to solve a complicated PDE for a function $h(x, y,t)$ with mixed derivatives and also a term involving $(h_y)^2$.

Just to try and get my basics right first, suppose we consider ( for $h(y,t)$ ) in a terminal value problem ( means we're solving backwards in time ),

$$h_t = (h_y)^2$$

For explicit Euler method we have ( backward differences in time and central differences in space ), $$\frac{h_n^m-h_n^{m-1}}{\Delta t} = \frac{(h_{n+1}^m)^2-2h_{n+1}^mh_{n-1}^m+(h_{n-1}^m)^2}{4\Delta y^2}\\ h_n^{m-1}=h_n^m - \frac{\Delta t}{4\Delta y^2}[(h_{n+1}^m)^2-2h_{n+1}^mh_{n-1}^m+(h_{n-1}^m)^2]$$ For this do I use something like von Neumann analysis?

For implicit Euler method we have ( forward differences in time and central differences in space ), $$\frac{h_n^{m+1}-h_n^{m}}{\Delta t} = \frac{(h_{n+1}^m)^2-2h_{n+1}^mh_{n-1}^m+(h_{n-1}^m)^2}{4\Delta y^2}\\ h_n^{m+1}=h_n^m + \frac{\Delta t}{4\Delta y^2}[(h_{n+1}^m)^2-2h_{n+1}^mh_{n-1}^m+(h_{n-1}^m)^2]$$ In this case we need to solve a system of nonlinear equations using Newton Rhapson or something like that.

The actual PDE i'm trying to solve is something more complicated which looks like this,

$$h_t + \left(\kappa(\theta-y)-\rho\mu\xi\right)h_y + \frac{1}{2}\xi^2yh_{yy} + \frac{1}{2}\gamma(1-\rho^2)\xi^2y h_y^2 + \xi^2y\frac{h_yf_y}{f} \\+ \frac{1}{2}yS^2h_{SS} + \rho\xi y S h_{Sy} =0\\h(S,y,T)= g(S,y)$$ where $\kappa, \theta,\xi, \gamma,\rho$ are constants, and $f$ is another given function.

I would very much appreciate if someone help with pointing me to good references/papers...

Thanks!

• Do you have spatial boundary conditions for $y$ and $S$? If either $|\rho|\geq 1$, $y\leq 0$, or $S=0$ then the linearized equation is not parabolic and this will be very difficult to accurately solve. – Aaron Oct 18 '15 at 1:41
• My $g(S,y)=g(S)$, which I made equal to an approximation to the heavy side function using $h(S,y,T)=g(S)=1/(1+e^{-\alpha(S-0.5)})$. – Danny Oct 19 '15 at 3:45
• I made some progress using an implicit finite difference and then sticking the nonlinear system into a Newton-Krylov solver available in scipy. It managed to evolve the solution back from T up to 17k/65k time steps, after which the solver got stuck because the jacobian matrix got all zeros and it wasn't able to descend anymore. – Danny Oct 19 '15 at 3:49
• @Aaron can you point me to references which describe how to linearize the squared derivative term? Will very much appreciate it – Danny Oct 19 '15 at 3:51
• You could read about the Fr\'{e}chet derivative in a decent textbook on functional analysis. This is just a Taylor expansion for non-linear operators on function spaces. If you wish to approximate the PDE about a function $\tilde{h}(S,y,t)$ then the non-linear term could be considered to be $(\tilde{h}_y+u_y)^2\approx \tilde{h}_y+2\tilde{h}_yu_y$. Since the other terms are linear in $h$, this would simply amount to replacing $h_y^2$ with $2\tilde{h}_yu_y$ (where $u_y$ is the variable you wish to solve for). I still have to ask what is the domain for $y$ and $S$ that you want to solve on? – Aaron Oct 20 '15 at 16:51