FEATool Multiphysics  v1.10
Finite Element Analysis Toolbox
Black-Scholes Equation

This is an example showing how to define a custom PDE equation model in the FEATool Multiphysics. In this case the Black-Scholes model equation, which is used in financial analytics to model derivatives and options pricing. The non-linear partial differential equation to be solved reads

\[ \frac{\partial u}{\partial t} - \frac{\partial}{\partial x}(\frac{1}{2}\frac{\partial u}{\partial x}) - \frac{\partial u}{\partial x} = -u + ( (x-t)^5 - 10(x-t)^4 - 10(x-t)^3) \]

with boundary conditions u(0,t) = -t^5 and u(x,t) = (x-t)^5 on the left and right sides of the domain, respectively. The problem is time-dependent with initial condition is u(x,t=0) = x^5. For this problem an exact analytical solution exists, u(x,t) = (x-t)^5, which can be used to verify the computed solution. See the Custom Equation section for details on the equation syntax.


This model is available as an automated tutorial by selecting Model Examples and Tutorials... > Classic PDE > Black-Scholes Equation from the File menu. Or alternatively, follow the step-by-step instructions below.

  1. To start a new model click the New Model toolbar button, or select New Model... from the File menu.
  2. Select the 1D radio button.
  3. Select the Custom Equation physics mode from the Select Physics drop-down menu.

  4. Press OK to finish the physics mode selection.
  5. Press the Create line Toolbar button.
  6. Press OK to finish and close the dialog box.

  7. Switch to Grid mode by clicking on the corresponding Mode Toolbar button.
  8. Press the Generate button to call the grid generation algorithm.

  9. Switch to Equation mode by clicking on the corresponding Mode Toolbar button.
  10. Press the edit button.
  11. Enter u' - 1/2*ux_x - ux_t + u_t = (x-t)^5 - 10*(x-t)^4 - 10*(x-t)^3 into the Equation for u edit field.

  12. Press OK to finish and close the dialog box.
  13. Enter x^5 into the Initial condition for u edit field.

  14. Press OK to finish the equation and subdomain settings specification.
  15. Switch to Boundary mode by clicking on the corresponding Mode Toolbar button.
  16. Enter -t^5 into the Dirichlet/Neumann coefficient edit field.

  17. Select 2 in the Boundaries list box.
  18. Enter (1-t)^5 into the Dirichlet/Neumann coefficient edit field.

  19. Press OK to finish the boundary condition specification.
  20. Switch to Solve mode by clicking on the corresponding Mode Toolbar button.
  21. Press the Settings Toolbar button.
  22. Select Time-Dependent from the Solution and solver type drop-down menu.

  23. Press the Solve button.

  24. To visualize the difference between the computed and exact analytical solution, open the Plot Options and postprocessing settings dialog box and enter the expression u - (x-t)^5 in the Surface Plot expression edit field.

  25. As the magnitude of the error too small to see clearly for the default axis scaling, manually zoom by selecting the Axis/Grid Settings... in the Options menu.
  26. Select the Axis equal radio button.
  27. Clear the Bounding box check box.
  28. Select the Axis limits radio button.
  29. Enter 0 1 -5e-3 5e-3 into the Manual axis settings [xmin xmax ymin ymax] edit field.

  30. Press OK to finish and close the dialog box.


One can now clearly see that the maximum error has a magnitude of 2e-3 around x = 0.25 which can improved further by using a finer grid and smaller time steps.

The black-scholes equation classic pde model has now been completed and can be saved as a binary (.fea) model file, or exported as a programmable MATLAB m-script text file, or GUI script (.fes) file.