a. The finite differences approach.

- Description of the method.

    The description of the finite differences approch could be found on this link.

- The reduction of the Black-Scholes equation.

    The Black-Scholes equation is given as defined previously :
    The term before the second derivative is positive, formally it lead to a negative diffusion. A steady perturbation, exp(-k*S), will be amplified by the equation an numerical instabilities will result.
    The negative coefficient behind V is also a cause of divergence, lets consider a uniform perturbation, exp(a*t). This equation gives that a=r>0. The perturbation will also be amplified and numerical instabilities could raise. To avoid these problem, the resolution of a modified equation is a elegant way to avoid these problems. The following changements, a new time varable

with t contained in [0, Ts2/2], a new space varable

with x between minus infinity and 0.47 if E=10 and

a new function u(x,t)
V(S,t)=
where
k1=
lead to a basic diffusion equation :

- Numerical application.

The transform equation has been solved with an implicit finite difference, 2 level scheme which is quite successful for such diffusion equations. the results are running in the following applet.
 

align=center width=800 height=280> 
 

Comparison with resolution using a Monte Carlo method is available on this link.
 

The resulting code, implemented in jbone is given below :
     //=======================================================================
      //  FD  - Black Scholes equations
      //=======================================================================
      // lets make some "short cut"
      BandMatrix a = new BandMatrix(3, f.length);
      BandMatrix b = new BandMatrix(3, f.length);
      double[]   c = new double[f.length];
      double[] cst = new double[f.length];

      // resolution of the diff eq as a tranformation of the BS eq

 // definition of some constants
        double step=jbone.step;
        double sigma = 0.45;
        double Dzero = 0.;
        double Te = 0.25;
        double re = 0.1;
        double Eu = 10.0;
        double kone = 2*re/(sigma*sigma);
        double ktwo = 2*(re-Dzero)/(sigma*sigma);
        double h=dx[0];
        double tempsAv=timeStep*(double)(step);
        double tauAv=(tempsAv)/2*sigma*sigma;
        double ix=0.;
        double[] dix = new double[f.length];
        double dtau=(timeStep)/2*sigma*sigma;
        double bdroit =5.75;
        double bgauche = 0.;

 if (step ==1) {bdroit=f[n-1];
 bgauche=f[0];}
 
      for (int i=1;i<=n; i++) {
        dix[i]=(Math.log(h*(double)(i+1.)/Eu))-(Math.log(h*(double)(i)/Eu));
      }
        dix[0]=(Math.log(h*(1.)/Eu))+10.;
 
      // Resolution of a basic diffusion equation

      for (int i=1; i<n; i++) {
        a.setL(i,  -2*dtau/dix[i-1]/(dix[i]+dix[i-1]));             //Matrix elements
        a.setD(i,   1.+2*dtau/dix[i]/dix[i-1]  );
        a.setR(i,  -2*dtau/dix[i]/(dix[i]+dix[i-1]));
        b.setL(i,   0.      );                                        //Right hand side
        b.setD(i,   1.      );
        b.setR(i,   0.      );
      }
       // boundary conditions
        a.setR(0, 0. );             //Matrix elements
        a.setD(0, 1. );
        a.setL(0, 0. );
        b.setL(0, 0. );             //Right hand side
        b.setD(0, 1. );
        b.setR(0, 0. );
 
        b.setL(n, 0. );             //Right hand side
        b.setD(n, 1. );
        b.setR(n, 0. );
        a.setD(n, 1. );
        a.setL(n, 0. );

      c=b.dot(g);                                    //Right hand side
            c[0]=bgauche/(Eu*Math.exp(-0.5*(ktwo-1.)*0.-(.25*(ktwo-1.)*(ktwo-1.)*0.+kone)*tauAv));                      // with periodicity
      c[n]=bdroit/(Eu*Math.exp(-0.5*(ktwo-1.)*0.47-(.25*(ktwo-1.)*(ktwo-1.)*0.47+kone)*tauAv));
      gp=a.solve3(c);

      //conversion of u -> V (return to the Sspace to get V)
 
      for (int i=1; i<=n; i++) {
      ix=Math.log((double)(i)*h/Eu);
      fp[i]=gp[i]*(Eu*
         Math.exp(-0.5*(ktwo-1.)*ix-(.25*(ktwo-1.)*(ktwo-1.)*ix+kone)*tauAv));
      }
 
      fp[0]=f[0];
 
 

      isDefined = true;