//This code assumes relatively large diffusion coefficient (I used 0.3) //if you make it small, than well-known oscillatory effect arise //due to central difference approximartion of the first //order derivative, so keep this in mind and if you want small vall, than //change approximation of the first order derivative //or use "particle methods" //also you need "dumka.c" to run this //depending on your compiler you can attach it in different manner //the code suppose tobeself-explanatory,rather than very efficient, but if //you can sqeese more than 20% efficiency out of it-let me know please. #include "stdafx.h"//commen this line out if you are not unver Miceosoft Visual Studio #include #include #ifndef max #define max(a,b) #endif FILE *payoffFunction; FILE *solution; FILE *solution_exact; double t_old; int nPointsX=300; double x_max=1.5; double x_min=0.0; double intervalX=x_max-x_min; long counter_call_back=0; double mean_Reverting=0.05e0; //double r0=0.09e0;//0.09e0; double strike=0.025e0; double rate_of_meanreversion=2.0; double pi=acos(-1.0e0); //put instanteneouse square of vol here, // for example use 0.3*0.3*x // for SDE: // dx = a(x,t)dt + 0.3 *sqrt(x)*dW // double sigma_sq(double t, double x){ double sigma=0.3e0; return sigma*sigma; } // //right hand side goes here dV/dt = d(x0-x)*v/dx + d^2 sigma^2*v/d x^2 // void rhs(long neqn, double* y, double t,double *result){ double y_Lft, y_Rgt; double x; double dx=intervalX/(double) (nPointsX+1); double dx_sq_inv=1.0e0/(dx*dx); double dx_inv=1.0e0/dx; double x_Rgt, x_Left; //calculation of the right hand side by method of lines //also read comments above about sigma_sq-it got to be large //if you sigma is small change approximation of the first derivative //to avoid well-known oscillatory effect for( int i=0; i < neqn; i++ ) { x=x_min+dx*(double)i; x_Rgt=x+dx; x_Left=x-dx; if( i==0 ) { y_Lft=2.0e0*y[i]-y[i+1];//null gamma boundary } else { y_Lft=y[i-1]; } if( i==neqn-1 ){ y_Rgt=2.0e0*y[i]-y[i-1];//zero second derivative boundary } else { y_Rgt = y[i+1]; } result[i]=rate_of_meanreversion*(mean_Reverting-x)*(y_Rgt- y_Lft)*(0.5e0*dx_inv) + 0.5e0*sigma_sq(t,x)*(y_Rgt - 2.0e0*y[i]+y_Lft)*dx_sq_inv - x*y[i]; } //this part is for output only. //If you want efficiency-comment it out if ( counter_call_back%1 == 0 && t_old < t ) { t_old=t; double tau=1.0e0-t; fprintf(solution, "%f",tau ); for ( int m=0; m < neqn; m++ ) { x=x_min+dx*(double)m; fprintf( solution, " %f",y[m]); }; fprintf(solution,"\n"); //below is output of exact solution fprintf(solution_exact, "%f", tau ); for ( int i=0; i