step-78.h

Go to the documentation of this file.

) const
 *   {
 * #ifdef MMS
 *     return -Utilities::fixed_power<2, double>(this->get_time()) + 6;
 * #else
 *     return 0.;
 * #endif
 *   }
 * 
 * 
 * 
 * @endcode
 * 
 * Then, we handle the right boundary condition.
 * 
 * @code
 *   template <int dim>
 *   class RightBoundaryValues : public Function<dim>
 *   {
 *   public:
 *     RightBoundaryValues(const double strike_price, const double interest_rate);
 * 
 *     virtual double value(const Point<dim> & p,
 *                          const unsigned int component = 0) const override;
 * 
 *   private:
 *     const double strike_price;
 *     const double interest_rate;
 *   };
 * 
 * 
 *   template <int dim>
 *   RightBoundaryValues<dim>::RightBoundaryValues(const double strike_price,
 *                                                 const double interest_rate)
 *     : strike_price(strike_price)
 *     , interest_rate(interest_rate)
 *   {}
 * 
 * 
 *   template <int dim>
 *   double RightBoundaryValues<dim>::value(const Point<dim> & p,
 *                                          const unsigned int component) const
 *   {
 * #ifdef MMS
 *     return -Utilities::fixed_power<2, double>(p(component)) -
 *            Utilities::fixed_power<2, double>(this->get_time()) + 6;
 * #else
 *     return (p(component) - strike_price) *
 *            exp((-interest_rate) * (this->get_time()));
 * #endif
 *   }
 * 
 * 
 * 
 * @endcode
 * 
 * Finally, we handle the right hand side.
 * 
 * @code
 *   template <int dim>
 *   class RightHandSide : public Function<dim>
 *   {
 *   public:
 *     RightHandSide(const double asset_volatility, const double interest_rate);
 * 
 *     virtual double value(const Point<dim> & p,
 *                          const unsigned int component = 0) const override;
 * 
 *   private:
 *     const double asset_volatility;
 *     const double interest_rate;
 *   };
 * 
 * 
 *   template <int dim>
 *   RightHandSide<dim>::RightHandSide(const double asset_volatility,
 *                                     const double interest_rate)
 *     : asset_volatility(asset_volatility)
 *     , interest_rate(interest_rate)
 *   {}
 * 
 * 
 *   template <int dim>
 *   double RightHandSide<dim>::value(const Point<dim> & p,
 *                                    const unsigned int component) const
 *   {
 * #ifdef MMS
 *     return 2 * (this->get_time()) -
 *            Utilities::fixed_power<2, double>(asset_volatility * p(component)) -
 *            2 * interest_rate * Utilities::fixed_power<2, double>(p(component)) -
 *            interest_rate *
 *              (-Utilities::fixed_power<2, double>(p(component)) -
 *               Utilities::fixed_power<2, double>(this->get_time()) + 6);
 * #else
 *     (void)p;
 *     (void)component;
 *     return 0.0;
 * #endif
 *   }
 * 
 * 
 * 
 * @endcode
 * 
 * 
 * <a name="ThecodeBlackScholescodeClass"></a> 
 * <h3>The <code>BlackScholes</code> Class</h3>
 * 
 
 * 
 * The next piece is the declaration of the main class of this program. This
 * is very similar to the @ref step_26 "step-26" tutorial, with some modifications. New
 * matrices had to be added to calculate the A and B matrices, as well as the
 * @f$V_{diff}@f$ vector mentioned in the introduction. We also define the
 * parameters used in the problem.
 *   
 
 * 
 * - <code>maximum_stock_price</code>: The imposed upper bound on the spatial
 * domain. This is the maximum allowed stock price.
 * - <code>maturity_time</code>: The upper bound on the time domain. This is
 * when the option expires.\n
 * - <code>asset_volatility</code>: The volatility of the stock price.\n
 * - <code>interest_rate</code>: The risk free interest rate.\n
 * - <code>strike_price</code>: The agreed upon price that the buyer will
 * have the option of purchasing  the stocks at the expiration time.
 *   
 
 * 
 * Some slight differences between this program and @ref step_26 "step-26" are the creation
 * of the <code>a_matrix</code> and the <code>b_matrix</code>, which is
 * described in the introduction. We then also need to store the current time,
 * the size of the time step, and the number of the current time step.
 * Next, we will store the output into a <code>DataOutStack</code>
 * variable because we will be layering the solution at each time on top of
 * one another to create the solution manifold. Then, we have a variable that
 * stores the current cycle and number of cycles that we will run when
 * calculating the solution. The cycle is one full solution calculation given
 * a mesh. We refine the mesh once in between each cycle to exhibit the
 * convergence properties of our program. Finally, we store the convergence
 * data into a convergence table.
 *   
 
 * 
 * As far as member functions are concerned, we have a function that
 * calculates the convergence information for each cycle, called
 * <code>process_solution</code>. This is just like what is done in @ref step_7 "step-7".
 * 
 * @code
 *   template <int dim>
 *   class BlackScholes
 *   {
 *   public:
 *     BlackScholes();
 * 
 *     void run();
 * 
 *   private:
 *     void setup_system();
 *     void solve_time_step();
 *     void refine_grid();
 *     void process_solution();
 *     void add_results_for_output();
 *     void write_convergence_table();
 * 
 *     const double maximum_stock_price;
 *     const double maturity_time;
 *     const double asset_volatility;
 *     const double interest_rate;
 *     const double strike_price;
 * 
 *     Triangulation<dim> triangulation;
 *     FE_Q<dim>          fe;
 *     DoFHandler<dim>    dof_handler;
 * 
 *     AffineConstraints<double> constraints;
 * 
 *     SparsityPattern      sparsity_pattern;
 *     SparseMatrix<double> mass_matrix;
 *     SparseMatrix<double> laplace_matrix;
 *     SparseMatrix<double> a_matrix;
 *     SparseMatrix<double> b_matrix;
 *     SparseMatrix<double> system_matrix;
 * 
 *     Vector<double> solution;
 *     Vector<double> system_rhs;
 * 
 *     double time;
 *     double time_step;
 * 
 *     const double       theta;
 *     const unsigned int n_cycles;
 *     const unsigned int n_time_steps;
 * 
 *     DataOutStack<dim>        data_out_stack;
 *     std::vector<std::string> solution_names;
 * 
 *     ConvergenceTable convergence_table;
 *   };
 * 
 * @endcode
 * 
 * 
 * <a name="ThecodeBlackScholescodeImplementation"></a> 
 * <h3>The <code>BlackScholes</code> Implementation</h3>
 * 
 
 * 
 * Now, we get to the implementation of the main class. We will set the values
 * for the various parameters used in the problem. These were chosen because
 * they are fairly normal values for these parameters. Although the stock
 * price has no upper bound in reality (it is in fact infinite), we impose
 * an upper bound that is twice the strike price. This is a somewhat arbitrary
 * choice to be twice the strike price, but it is large enough to see the
 * interesting parts of the solution.
 * 
 * @code
 *   template <int dim>
 *   BlackScholes<dim>::BlackScholes()
 *     : maximum_stock_price(1.)
 *     , maturity_time(1.)
 *     , asset_volatility(.2)
 *     , interest_rate(0.05)
 *     , strike_price(0.5)
 *     , fe(1)
 *     , dof_handler(triangulation)
 *     , time(0.0)
 *     , theta(0.5)
 *     , n_cycles(4)
 *     , n_time_steps(5000)
 *   {
 *     Assert(dim == 1, ExcNotImplemented());
 *   }
 * 
 * @endcode
 * 
 * 
 * <a name="codeBlackScholessetup_systemcode"></a> 
 * <h4><code>BlackScholes::setup_system</code></h4>
 * 
 
 * 
 * The next function sets up the DoFHandler object, computes
 * the constraints, and sets the linear algebra objects to their correct
 * sizes. We also compute the mass matrix here by calling a function from the
 * library. We will compute the other 3 matrices next, because these need to
 * be computed 'by hand'.
 *   
 
 * 
 * Note, that the time step is initialized here because the maturity time was
 * needed to compute the time step.
 * 
 * @code
 *   template <int dim>
 *   void BlackScholes<dim>::setup_system()
 *   {
 *     dof_handler.distribute_dofs(fe);
 * 
 *     time_step = maturity_time / n_time_steps;
 * 
 *     constraints.clear();
 *     DoFTools::make_hanging_node_constraints(dof_handler, constraints);
 *     constraints.close();
 *     DynamicSparsityPattern dsp(dof_handler.n_dofs());
 *     DoFTools::make_sparsity_pattern(dof_handler,
 *                                     dsp,
 *                                     constraints,
 *                                     /*keep_constrained_dofs = */ true);
 *     sparsity_pattern.copy_from(dsp);
 * 
 *     mass_matrix.reinit(sparsity_pattern);
 *     laplace_matrix.reinit(sparsity_pattern);
 *     a_matrix.reinit(sparsity_pattern);
 *     b_matrix.reinit(sparsity_pattern);
 *     system_matrix.reinit(sparsity_pattern);
 * 
 *     MatrixCreator::create_mass_matrix(dof_handler,
 *                                       QGauss<dim>(fe.degree + 1),
 *                                       mass_matrix);
 * 
 * @endcode
 * 
 * Below is the code to create the Laplace matrix with non-constant
 * coefficients. This corresponds to the matrix D in the introduction. This
 * non-constant coefficient is represented in the
 * <code>current_coefficient</code> variable.
 * 
 * @code
 *     const unsigned int dofs_per_cell = fe.dofs_per_cell;
 *     FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
 *     QGauss<dim>        quadrature_formula(fe.degree + 1);
 *     FEValues<dim>      fe_values(fe,
 *                             quadrature_formula,
 *                             update_values | update_gradients |
 *                               update_quadrature_points | update_JxW_values);
 *     std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
 *     for (const auto &cell : dof_handler.active_cell_iterators())
 *       {
 *         cell_matrix = 0.;
 *         fe_values.reinit(cell);
 *         for (const unsigned int q_index : fe_values.quadrature_point_indices())
 *           {
 *             const double current_coefficient =
 *               fe_values.quadrature_point(q_index).square();
 *             for (const unsigned int i : fe_values.dof_indices())
 *               {
 *                 for (const unsigned int j : fe_values.dof_indices())
 *                   cell_matrix(i, j) +=
 *                     (current_coefficient *              // (x_q)^2
 *                      fe_values.shape_grad(i, q_index) * // grad phi_i(x_q)
 *                      fe_values.shape_grad(j, q_index) * // grad phi_j(x_q)
 *                      fe_values.JxW(q_index));           // dx
 *               }
 *           }
 *         cell->get_dof_indices(local_dof_indices);
 *         for (const unsigned int i : fe_values.dof_indices())
 *           {
 *             for (const unsigned int j : fe_values.dof_indices())
 *               laplace_matrix.add(local_dof_indices[i],
 *                                  local_dof_indices[j],
 *                                  cell_matrix(i, j));
 *           }
 *       }
 * 
 * @endcode
 * 
 * Now we will create the A matrix. Below is the code to create the matrix A
 * as discussed in the introduction. The non constant coefficient is again
 * represented in  the <code>current_coefficient</code> variable.
 * 
 * @code
 *     for (const auto &cell : dof_handler.active_cell_iterators())
 *       {
 *         cell_matrix = 0.;
 *         fe_values.reinit(cell);
 *         for (const unsigned int q_index : fe_values.quadrature_point_indices())
 *           {
 *             const Tensor<1, dim> current_coefficient =
 *               fe_values.quadrature_point(q_index);
 *             for (const unsigned int i : fe_values.dof_indices())
 *               {
 *                 for (const unsigned int j : fe_values.dof_indices())
 *                   {
 *                     cell_matrix(i, j) +=
 *                       (current_coefficient *               // x_q
 *                        fe_values.shape_grad(i, q_index) *  // grad phi_i(x_q)
 *                        fe_values.shape_value(j, q_index) * // phi_j(x_q)
 *                        fe_values.JxW(q_index));            // dx
 *                   }
 *               }
 *           }
 *         cell->get_dof_indices(local_dof_indices);
 *         for (const unsigned int i : fe_values.dof_indices())
 *           {
 *             for (const unsigned int j : fe_values.dof_indices())
 *               a_matrix.add(local_dof_indices[i],
 *                            local_dof_indices[j],
 *                            cell_matrix(i, j));
 *           }
 *       }
 * 
 * @endcode
 * 
 * Finally we will create the matrix B. Below is the code to create the
 * matrix B as discussed in the introduction. The non constant coefficient
 * is again represented in the <code>current_coefficient</code> variable.
 * 
 * @code
 *     for (const auto &cell : dof_handler.active_cell_iterators())
 *       {
 *         cell_matrix = 0.;
 *         fe_values.reinit(cell);
 *         for (const unsigned int q_index : fe_values.quadrature_point_indices())
 *           {
 *             const Tensor<1, dim> current_coefficient =
 *               fe_values.quadrature_point(q_index);
 *             for (const unsigned int i : fe_values.dof_indices())
 *               {
 *                 for (const unsigned int j : fe_values.dof_indices())
 *                   cell_matrix(i, j) +=
 *                     (current_coefficient *               // x_q
 *                      fe_values.shape_value(i, q_index) * // phi_i(x_q)
 *                      fe_values.shape_grad(j, q_index) *  // grad phi_j(x_q)
 *                      fe_values.JxW(q_index));            // dx
 *               }
 *           }
 *         cell->get_dof_indices(local_dof_indices);
 *         for (const unsigned int i : fe_values.dof_indices())
 *           {
 *             for (const unsigned int j : fe_values.dof_indices())
 *               b_matrix.add(local_dof_indices[i],
 *                            local_dof_indices[j],
 *                            cell_matrix(i, j));
 *           }
 *       }
 * 
 *     solution.reinit(dof_handler.n_dofs());
 *     system_rhs.reinit(dof_handler.n_dofs());
 *   }
 * 
 * @endcode
 * 
 * 
 * <a name="codeBlackScholessolve_time_stepcode"></a> 
 * <h4><code>BlackScholes::solve_time_step</code></h4>
 * 
 
 * 
 * The next function is the one that solves the actual linear system for a
 * single time step. The only interesting thing here is that the matrices
 * we have built are symmetric positive definite, so we can use the
 * conjugate gradient method.
 * 
 * @code
 *   template <int dim>
 *   void BlackScholes<dim>::solve_time_step()
 *   {
 *     SolverControl                          solver_control(1000, 1e-12);
 *     SolverCG<Vector<double>>               cg(solver_control);
 *     PreconditionSSOR<SparseMatrix<double>> preconditioner;
 *     preconditioner.initialize(system_matrix, 1.0);
 *     cg.solve(system_matrix, solution, system_rhs, preconditioner);
 *     constraints.distribute(solution);
 *   }
 * 
 * @endcode
 * 
 * 
 * <a name="codeBlackScholesadd_results_for_outputcode"></a> 
 * <h4><code>BlackScholes::add_results_for_output</code></h4>
 * 
 
 * 
 * This is simply the function to stitch the solution pieces together. For
 * this, we create a new layer at each time, and then add the solution vector
 * for that timestep. The function then stitches this together with the old
 * solutions using 'build_patches'.
 * 
 * @code
 *   template <int dim>
 *   void BlackScholes<dim>::add_results_for_output()
 *   {
 *     data_out_stack.new_parameter_value(time, time_step);
 *     data_out_stack.attach_dof_handler(dof_handler);
 *     data_out_stack.add_data_vector(solution, solution_names);
 *     data_out_stack.build_patches(2);
 *     data_out_stack.finish_parameter_value();
 *   }
 * 
 * @endcode
 * 
 * 
 * <a name="codeBlackScholesrefine_gridcode"></a> 
 * <h4><code>BlackScholes::refine_grid</code></h4>
 * 
 
 * 
 * It is somewhat unnecessary to have a function for the global refinement
 * that we do. The reason for the function is to allow for the possibility of
 * an adaptive refinement later.
 * 
 * @code
 *   template <int dim>
 *   void BlackScholes<dim>::refine_grid()
 *   {
 *     triangulation.refine_global(1);
 *   }
 * 
 * @endcode
 * 
 * 
 * <a name="codeBlackScholesprocess_solutioncode"></a> 
 * <h4><code>BlackScholes::process_solution</code></h4>
 * 
 
 * 
 * This is where we calculate the convergence and error data to evaluate the
 * effectiveness of the program. Here, we calculate the @f$L^2@f$, @f$H^1@f$ and
 * @f$L^{\infty}@f$ norms.
 * 
 * @code
 *   template <int dim>
 *   void BlackScholes<dim>::process_solution()
 *   {
 *     Solution<dim> sol(maturity_time);
 *     sol.set_time(time);
 *     Vector<float> difference_per_cell(triangulation.n_active_cells());
 *     VectorTools::integrate_difference(dof_handler,
 *                                       solution,
 *                                       sol,
 *                                       difference_per_cell,
 *                                       QGauss<dim>(fe.degree + 1),
 *                                       VectorTools::L2_norm);
 *     const double L2_error =
 *       VectorTools::compute_global_error(triangulation,
 *                                         difference_per_cell,
 *                                         VectorTools::L2_norm);
 *     VectorTools::integrate_difference(dof_handler,
 *                                       solution,
 *                                       sol,
 *                                       difference_per_cell,
 *                                       QGauss<dim>(fe.degree + 1),
 *                                       VectorTools::H1_seminorm);
 *     const double H1_error =
 *       VectorTools::compute_global_error(triangulation,
 *                                         difference_per_cell,
 *                                         VectorTools::H1_seminorm);
 *     const QTrapezoid<1>  q_trapezoid;
 *     const QIterated<dim> q_iterated(q_trapezoid, fe.degree * 2 + 1);
 *     VectorTools::integrate_difference(dof_handler,
 *                                       solution,
 *                                       sol,
 *                                       difference_per_cell,
 *                                       q_iterated,
 *                                       VectorTools::Linfty_norm);
 *     const double Linfty_error =
 *       VectorTools::compute_global_error(triangulation,
 *                                         difference_per_cell,
 *                                         VectorTools::Linfty_norm);
 *     const unsigned int n_active_cells = triangulation.n_active_cells();
 *     const unsigned int n_dofs         = dof_handler.n_dofs();
 *     convergence_table.add_value("cells", n_active_cells);
 *     convergence_table.add_value("dofs", n_dofs);
 *     convergence_table.add_value("L2", L2_error);
 *     convergence_table.add_value("H1", H1_error);
 *     convergence_table.add_value("Linfty", Linfty_error);
 *   }
 * 
 * @endcode
 * 
 * 
 * <a name="codeBlackScholeswrite_convergence_tablecode"></a> 
 * <h4><code>BlackScholes::write_convergence_table</code> </h4>
 * 
 
 * 
 * This next part is building the convergence and error tables. By this, we
 * need to set the settings for how to output the data that was calculated
 * during <code>BlackScholes::process_solution</code>. First, we will create
 * the headings and set up the cells properly. During this, we will also
 * prescribe the precision of our results. Then we will write the calculated
 * errors based on the @f$L^2@f$, @f$H^1@f$, and @f$L^{\infty}@f$ norms to the console and
 * to the error LaTeX file.
 * 
 * @code
 *   template <int dim>
 *   void BlackScholes<dim>::write_convergence_table()
 *   {
 *     convergence_table.set_precision("L2", 3);
 *     convergence_table.set_precision("H1", 3);
 *     convergence_table.set_precision("Linfty", 3);
 *     convergence_table.set_scientific("L2", true);
 *     convergence_table.set_scientific("H1", true);
 *     convergence_table.set_scientific("Linfty", true);
 *     convergence_table.set_tex_caption("cells", "\\# cells");
 *     convergence_table.set_tex_caption("dofs", "\\# dofs");
 *     convergence_table.set_tex_caption("L2", "@fL^2@f-error");
 *     convergence_table.set_tex_caption("H1", "@fH^1@f-error");
 *     convergence_table.set_tex_caption("Linfty", "@fL^\\infty@f-error");
 *     convergence_table.set_tex_format("cells", "r");
 *     convergence_table.set_tex_format("dofs", "r");
 *     std::cout << std::endl;
 *     convergence_table.write_text(std::cout);
 *     std::string error_filename = "error";
 *     error_filename += "-global";
 *     error_filename += ".tex";
 *     std::ofstream error_table_file(error_filename);
 *     convergence_table.write_tex(error_table_file);
 * 
 * @endcode
 * 
 * Next, we will make the convergence table. We will again write this to
 * the console and to the convergence LaTeX file.
 * 
 * @code
 *     convergence_table.add_column_to_supercolumn("cells", "n cells");
 *     std::vector<std::string> new_order;
 *     new_order.emplace_back("n cells");
 *     new_order.emplace_back("H1");
 *     new_order.emplace_back("L2");
 *     convergence_table.set_column_order(new_order);
 *     convergence_table.evaluate_convergence_rates(
 *       "L2", ConvergenceTable::reduction_rate);
 *     convergence_table.evaluate_convergence_rates(
 *       "L2", ConvergenceTable::reduction_rate_log2);
 *     convergence_table.evaluate_convergence_rates(
 *       "H1", ConvergenceTable::reduction_rate);
 *     convergence_table.evaluate_convergence_rates(
 *       "H1", ConvergenceTable::reduction_rate_log2);
 *     std::cout << std::endl;
 *     convergence_table.write_text(std::cout);
 *     std::string conv_filename = "convergence";
 *     conv_filename += "-global";
 *     switch (fe.degree)
 *       {
 *         case 1:
 *           conv_filename += "-q1";
 *           break;
 *         case 2:
 *           conv_filename += "-q2";
 *           break;
 *         default:
 *           Assert(false, ExcNotImplemented());
 *       }
 *     conv_filename += ".tex";
 *     std::ofstream table_file(conv_filename);
 *     convergence_table.write_tex(table_file);
 *   }
 * 
 * @endcode
 * 
 * 
 * <a name="codeBlackScholesruncode"></a> 
 * <h4><code>BlackScholes::run</code></h4>
 * 
 
 * 
 * Now we get to the main driver of the program. This is where we do all the
 * work of looping through the time steps and calculating the solution vector
 * each time. Here at the top, we set the initial refinement value and then
 * create a mesh. Then we refine this mesh once. Next, we set up the
 * data_out_stack object to store our solution. Finally, we start a for loop
 * to loop through the cycles. This lets us recalculate a solution for each
 * successive mesh refinement. At the beginning of each iteration, we need to
 * reset the time and time step number. We introduce an if statement to
 * accomplish this because we don't want to do this on the first iteration.
 * 
 * @code
 *   template <int dim>
 *   void BlackScholes<dim>::run()
 *   {
 *     GridGenerator::hyper_cube(triangulation, 0.0, maximum_stock_price, true);
 *     triangulation.refine_global(0);
 * 
 *     solution_names.emplace_back("u");
 *     data_out_stack.declare_data_vector(solution_names,
 *                                        DataOutStack<dim>::dof_vector);
 * 
 *     Vector<double> vmult_result;
 *     Vector<double> forcing_terms;
 * 
 *     for (unsigned int cycle = 0; cycle < n_cycles; cycle++)
 *       {
 *         if (cycle != 0)
 *           {
 *             refine_grid();
 *             time = 0.0;
 *           }
 * 
 *         setup_system();
 * 
 *         std::cout << std::endl
 *                   << "===========================================" << std::endl
 *                   << "Cycle " << cycle << ':' << std::endl
 *                   << "Number of active cells: "
 *                   << triangulation.n_active_cells() << std::endl
 *                   << "Number of degrees of freedom: " << dof_handler.n_dofs()
 *                   << std::endl
 *                   << std::endl;
 * 
 *         VectorTools::interpolate(dof_handler,
 *                                  InitialConditions<dim>(strike_price),
 *                                  solution);
 * 
 *         if (cycle == (n_cycles - 1))
 *           {
 *             add_results_for_output();
 *           }
 * 
 * @endcode
 * 
 * Next, we run the main loop which runs until we exceed the maturity
 * time. We first compute the right hand side of the equation, which is
 * described in the introduction. Recall that it contains the term
 * @f$\left[-\frac{1}{4}k_n\sigma^2\mathbf{D}-k_nr\mathbf{M}+k_n\sigma^2
 * \mathbf{B}-k_nr\mathbf{A}+\mathbf{M}\right]V^{n-1}@f$. We put these
 * terms into the variable system_rhs, with the help of a temporary
 * vector:
 * 
 * @code
 *         vmult_result.reinit(dof_handler.n_dofs());
 *         forcing_terms.reinit(dof_handler.n_dofs());
 *         for (unsigned int timestep_number = 0; timestep_number < n_time_steps;
 *              ++timestep_number)
 *           {
 *             time += time_step;
 * 
 *             if (timestep_number % 1000 == 0)
 *               std::cout << "Time step " << timestep_number << " at t=" << time
 *                         << std::endl;
 * 
 *             mass_matrix.vmult(system_rhs, solution);
 * 
 *             laplace_matrix.vmult(vmult_result, solution);
 *             system_rhs.add(
 *               (-1) * (1 - theta) * time_step *
 *                 Utilities::fixed_power<2, double>(asset_volatility) * 0.5,
 *               vmult_result);
 *             mass_matrix.vmult(vmult_result, solution);
 * 
 *             system_rhs.add((-1) * (1 - theta) * time_step * interest_rate * 2,
 *                            vmult_result);
 * 
 *             a_matrix.vmult(vmult_result, solution);
 *             system_rhs.add((-1) * time_step * interest_rate, vmult_result);
 * 
 *             b_matrix.vmult(vmult_result, solution);
 *             system_rhs.add(
 *               (-1) * Utilities::fixed_power<2, double>(asset_volatility) *
 *                 time_step * 1,
 *               vmult_result);
 * 
 * @endcode
 * 
 * The second piece is to compute the contributions of the source
 * terms. This corresponds to the term @f$-k_n\left[\frac{1}{2}F^{n-1}
 * +\frac{1}{2}F^n\right]@f$. The following code calls
 * VectorTools::create_right_hand_side to compute the vectors @f$F@f$,
 * where we set the time of the right hand side (source) function
 * before we evaluate it. The result of this all ends up in the
 * forcing_terms variable:
 * 
 * @code
 *             RightHandSide<dim> rhs_function(asset_volatility, interest_rate);
 *             rhs_function.set_time(time);
 *             VectorTools::create_right_hand_side(dof_handler,
 *                                                 QGauss<dim>(fe.degree + 1),
 *                                                 rhs_function,
 *                                                 forcing_terms);
 *             forcing_terms *= time_step * theta;
 *             system_rhs -= forcing_terms;
 * 
 *             rhs_function.set_time(time - time_step);
 *             VectorTools::create_right_hand_side(dof_handler,
 *                                                 QGauss<dim>(fe.degree + 1),
 *                                                 rhs_function,
 *                                                 forcing_terms);
 *             forcing_terms *= time_step * (1 - theta);
 *             system_rhs -= forcing_terms;
 * 
 * @endcode
 * 
 * Next, we add the forcing terms to the ones that come from the
 * time stepping, and also build the matrix @f$\left[\mathbf{M}+
 * \frac{1}{4}k_n\sigma^2\mathbf{D}+k_nr\mathbf{M}\right]@f$ that we
 * have to invert in each time step. The final piece of these
 * operations is to eliminate hanging node constrained degrees of
 * freedom from the linear system:
 * 
 * @code
 *             system_matrix.copy_from(mass_matrix);
 *             system_matrix.add(
 *               (theta)*time_step *
 *                 Utilities::fixed_power<2, double>(asset_volatility) * 0.5,
 *               laplace_matrix);
 *             system_matrix.add((time_step)*interest_rate * theta * (1 + 1),
 *                               mass_matrix);
 * 
 *             constraints.condense(system_matrix, system_rhs);
 * 
 * @endcode
 * 
 * There is one more operation we need to do before we can solve it:
 * boundary values. To this end, we create a boundary value object,
 * set the proper time to the one of the current time step, and
 * evaluate it as we have done many times before. The result is used
 * to also set the correct boundary values in the linear system:
 * 
 * @code
 *             {
 *               RightBoundaryValues<dim> right_boundary_function(strike_price,
 *                                                                interest_rate);
 *               LeftBoundaryValues<dim>  left_boundary_function;
 *               right_boundary_function.set_time(time);
 *               left_boundary_function.set_time(time);
 *               std::map<types::global_dof_index, double> boundary_values;
 *               VectorTools::interpolate_boundary_values(dof_handler,
 *                                                        0,
 *                                                        left_boundary_function,
 *                                                        boundary_values);
 *               VectorTools::interpolate_boundary_values(dof_handler,
 *                                                        1,
 *                                                        right_boundary_function,
 *                                                        boundary_values);
 *               MatrixTools::apply_boundary_values(boundary_values,
 *                                                  system_matrix,
 *                                                  solution,
 *                                                  system_rhs);
 *             }
 * 
 * @endcode
 * 
 * With this out of the way, all we have to do is solve the system,
 * generate graphical data on the last cycle, and create the
 * convergence table data.
 * 
 * @code
 *             solve_time_step();
 * 
 *             if (cycle == (n_cycles - 1))
 *               {
 *                 add_results_for_output();
 *               }
 *           }
 * #ifdef MMS
 *         process_solution();
 * #endif
 *       }
 * 
 *     const std::string filename = "solution.vtk";
 *     std::ofstream     output(filename);
 *     data_out_stack.write_vtk(output);
 * 
 * #ifdef MMS
 *     write_convergence_table();
 * #endif
 *   }
 * 
 * } // namespace BlackScholesSolver
 * 
 * @endcode
 * 
 * 
 * <a name="ThecodemaincodeFunction"></a> 
 * <h3>The <code>main</code> Function</h3>
 * 
 
 * 
 * Having made it this far, there is, again, nothing much to discuss for the
 * main function of this program: it looks like all such functions since @ref step_6 "step-6".
 * 
 * @code
 * int main()
 * {
 *   try
 *     {
 *       using namespace BlackScholesSolver;
 * 
 *       BlackScholes<1> black_scholes_solver;
 *       black_scholes_solver.run();
 *     }
 *   catch (std::exception &exc)
 *     {
 *       std::cerr << std::endl
 *                 << std::endl
 *                 << "----------------------------------------------------"
 *                 << std::endl;
 *       std::cerr << "Exception on processing: " << std::endl
 *                 << exc.what() << std::endl
 *                 << "Aborting!" << std::endl
 *                 << "----------------------------------------------------"
 *                 << std::endl;
 *       return 1;
 *     }
 *   catch (...)
 *     {
 *       std::cerr << std::endl
 *                 << std::endl
 *                 << "----------------------------------------------------"
 *                 << std::endl;
 *       std::cerr << "Unknown exception!" << std::endl
 *                 << "Aborting!" << std::endl
 *                 << "----------------------------------------------------"
 *                 << std::endl;
 *       return 1;
 *     }
 *   return 0;
 * }
 * @endcode
<a name="Results"></a><h1>Results</h1>
 
 
 
Below is the output of the program:
@code
===========================================
Number of active cells: 1
Number of degrees of freedom: 2
 
Time step 0 at t=0.0002
[...]
 
Cycle 7:
Number of active cells: 128
Number of degrees of freedom: 129
 
Time step 0 at t=0.0002
Time step 1000 at t=0.2002
Time step 2000 at t=0.4002
Time step 3000 at t=0.6002
Time step 4000 at t=0.8002
 
cells dofs    L2        H1      Linfty
    1    2 1.667e-01 5.774e-01 2.222e-01
    2    3 3.906e-02 2.889e-01 5.380e-02
    4    5 9.679e-03 1.444e-01 1.357e-02
    8    9 2.405e-03 7.218e-02 3.419e-03
   16   17 5.967e-04 3.609e-02 8.597e-04
   32   33 1.457e-04 1.804e-02 2.155e-04
   64   65 3.307e-05 9.022e-03 5.388e-05
  128  129 5.016e-06 4.511e-03 1.342e-05
 
n cells         H1                  L2
      1 5.774e-01    -    - 1.667e-01    -    -
      2 2.889e-01 2.00 1.00 3.906e-02 4.27 2.09
      4 1.444e-01 2.00 1.00 9.679e-03 4.04 2.01
      8 7.218e-02 2.00 1.00 2.405e-03 4.02 2.01
     16 3.609e-02 2.00 1.00 5.967e-04 4.03 2.01
     32 1.804e-02 2.00 1.00 1.457e-04 4.10 2.03
     64 9.022e-03 2.00 1.00 3.307e-05 4.41 2.14
    128 4.511e-03 2.00 1.00 5.016e-06 6.59 2.72
@endcode
 
What is more interesting is the output of the convergence tables. They are
outputted into the console, as well into a LaTeX file. The convergence tables
are shown above. Here, you can see that the the solution has a convergence rate
of @f$\mathcal{O}(h)@f$ with respect to the @f$H^1@f$-norm, and the solution has a convergence rate
of @f$\mathcal{O}(h^2)@f$ with respect to the @f$L^2@f$-norm.
 
 
Below is the visualization of the solution.
 
<div style="text-align:center;">
  <img src="https://www.dealii.org/images/steps/developer/step-78.mms-solution.png"
       alt="Solution of the MMS problem.">
</div>
 *
 *
<a name="PlainProg"></a>
<h1> The plain program</h1>
@include "step-78.cc"
*/