b. The finite element 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 elements scheme with roof function and the integration is performe with a mid point approximation, 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 :

      //=======================================================================
      //  FEM  - Black Scholes equation
      //=======================================================================
      BandMatrix a = new BandMatrix(3, f.length);
      BandMatrix b = new BandMatrix(3, f.length);
      double[]   c = new double[f.length];
 

      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 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;

       // New coordinate frame
      for (int i=1;i<=n; i++) {
        dix[i]=(Math.log(dx[0]*(double)(i+1.)/Eu))-(Math.log(dx[0]*(double)(i)/Eu));
      }
        dix[0]=(Math.log(dx[0]/Eu))+10.;

 

       //Resolution of the diffusion equation
      for (int i=1; i<n; i++) {
        a.setR(i,                           dix[i]/4. -dtau/dix[i]* theta    );
        a.setD(i,1./4.*(dix[i]+dix[i-1])+(1./dix[i-1]+1./dix[i])*dtau* theta    );
        a.setL(i,                       dix[i-1]/4. -dtau/dix[i-1]* theta    );
        b.setR(i,                           dix[i]/4. -dtau/dix[i]*(theta-1) );
        b.setD(i,1./4.*(dix[i]+dix[i-1])+(1./dix[i-1]+1./dix[i])*dtau*(theta-1) );
        b.setL(i,                       dix[i-1]/4. -dtau/dix[i-1]*(theta-1) );
      }
 
 

         //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, 0. );
        b.setR(0, 0. );

 
         a.setD(n, 1. );
         a.setL(n, 0. );
         b.setL(n, 0. );             //Right hand side
         b.setD(n, 0. );
         b.setR(n, 0. );
 
      c=b.dot(g);                                 //Right hand side
      c[0]=0.;
 
      c[n]=5.75/(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 (from Xspace to Sspace)
 
      for (int i=1; i<=n; i++) {
      ix=Math.log((double)(i)*dx[0]/Eu);
      fp[i]=gp[i]*(Eu*
         Math.exp(-0.5*(ktwo-1.)*ix-(.25*(ktwo-1.)*(ktwo-1.)*ix+kone)*tauAv));
 
      }
      fp[0]=0.;
 
      isDefined = true;