diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/5.1IVP/IVP_test.md b/MATH5620_NumericalSolutionsOfDifferentialEquations/5.1IVP/IVP_test.md new file mode 100644 index 0000000..3e92002 --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/5.1IVP/IVP_test.md @@ -0,0 +1,81 @@ +--- +title: IVP Test +math: true +layout: default +--- + +{% include mathjax.html %} + + Table of Contents +# First Order IVP Solver + +**Routine Name:** firstOrderIVPSolver + +**Author:** Philip Nelson + +**Language:** C++ + +## Description + +This method is used to solve the simple first order initial value problem + +\\[u` = \lambda u\\] + +with \\(u(0)=P_0\\). This ordinary differential equations has the soltion + +\\[u(t)=\alpha e^{\lambda t}\\] + +## Input + +`firstOrderIVPSolver(const T& lambda, const T& alpha)` requires: + +* `const T& l` - \\(\lambda\\) +* `const T& a` - \\(\alpha\\) + +## Output + +This method returns a function that can be evaluated at any time \\(t\\) +to obtain the exact solution to the analytic equation +(including roundoff error). + +## Code +{% highlight c++ %} +template +auto firstOrderIVPSolver(const T& l, const T& a) +{ + return [=](const T& t) { return a * std::exp(l * t); }; +} +{% endhighlight %} + +## Example +{% highlight c++ %} +int main() +{ + + auto solveIVP = firstOrderIVPSolver(-1.5, 7.3); + + for(auto t = 0; t < 10; ++t) + { + std::cout << "t = " << t << " -> " << solveIVP(t) + << "\tt = " << t+10 << " -> " << solveIVP(t+10) << '\n'; + } + + return EXIT_SUCCESS; +} +{% endhighlight %} + +## Result +``` +t = 0 -> 7.3 t = 10 -> 2.23309e-06 +t = 1 -> 1.62885 t = 11 -> 4.98269e-07 +t = 2 -> 0.363446 t = 12 -> 1.11179e-07 +t = 3 -> 0.0810957 t = 13 -> 2.48074e-08 +t = 4 -> 0.0180949 t = 14 -> 5.53527e-09 +t = 5 -> 0.00403752 t = 15 -> 1.23509e-09 +t = 6 -> 0.000900892 t = 16 -> 2.75585e-10 +t = 7 -> 0.000201016 t = 17 -> 6.14913e-11 +t = 8 -> 4.48528e-05 t = 18 -> 1.37206e-11 +t = 9 -> 1.0008e-05 t = 19 -> 3.06147e-12 +``` + +**Last Modification date:** 3 April 2018 diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/5.1IVP/Makefile b/MATH5620_NumericalSolutionsOfDifferentialEquations/5.1IVP/Makefile new file mode 100644 index 0000000..83763f3 --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/5.1IVP/Makefile @@ -0,0 +1,12 @@ +OBJS = main.cpp firstOrderIVP.hpp + +CC = g++ +FLAGS = -std=c++17 +DEBUG_FLAGS = -Og -g3 -Wall -Wextra -Werror +RELEASE_FLAGS = -O3 + +release: $(OBJS) + $(CC) $(RELEASE_FLAGS) $(FLAGS) $(OBJS) -o ivp.out + +debug: $(OBJS) + $(CC) $(DEBUG_FLAGS) $(FLAGS) $(OBJS) -o debug.out diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/5.1IVP/firstOrderIVP.hpp b/MATH5620_NumericalSolutionsOfDifferentialEquations/5.1IVP/firstOrderIVP.hpp new file mode 100644 index 0000000..5d3922d --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/5.1IVP/firstOrderIVP.hpp @@ -0,0 +1,13 @@ +#ifndef FIRST_ORDER_IVP_HPP +#define FIRST_ORDER_IVP_HPP + +#include +#include + +template +auto firstOrderIVPSolver(const T& l, const T& a) +{ + return [=](const T& t) { return a * std::exp(l * t); }; +} + +#endif diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/5.1IVP/main.cpp b/MATH5620_NumericalSolutionsOfDifferentialEquations/5.1IVP/main.cpp new file mode 100644 index 0000000..7a6868c --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/5.1IVP/main.cpp @@ -0,0 +1,17 @@ +#include +#include +#include "firstOrderIVP.hpp" + +int main() +{ + + auto solveIVP = firstOrderIVPSolver(-1.5, 7.3); + + for(auto t = 0; t < 10; ++t) + { + std::cout << "t = " << t << " -> " << solveIVP(t) + << "\tt = " << t+10 << " -> " << solveIVP(t+10) << '\n'; + } + + return EXIT_SUCCESS; +} diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/HW/1/HW1-10.pdf b/MATH5620_NumericalSolutionsOfDifferentialEquations/HW/1/HW1-10.pdf deleted file mode 100644 index 901d341..0000000 Binary files a/MATH5620_NumericalSolutionsOfDifferentialEquations/HW/1/HW1-10.pdf and /dev/null differ diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/HW/1/eigen.html b/MATH5620_NumericalSolutionsOfDifferentialEquations/HW/1/eigen.html deleted file mode 100644 index 5478c35..0000000 --- a/MATH5620_NumericalSolutionsOfDifferentialEquations/HW/1/eigen.html +++ /dev/null @@ -1,420 +0,0 @@ - - - - - - - - - - - - - - -HW1.10 - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- - - - - - - - - - - - - - -
-

Compute the eigenvalues and eigenvectors for the following matrix using ye olde proverbally pencil and paper

-
from sympy import *
-A = Matrix([[3.4109, -0.2693, -1.0643],[1.5870, 1.5546, -5.5361],[0.2981, -0.2981, 1.2277]])
-I = Matrix([[1,0,0],[0,1,0],[0,0,1]])
-eVals = A.eigenvals()
-eVects = A.eigenvects()
-pprint(A)
-# A matrix
-
## [3.4109  -0.2693  -1.0643]
-## [                        ]
-## [1.587   1.5546   -5.5361]
-## [                        ]
-## [0.2981  -0.2981  1.2277 ]
-
for val in list(eVals.keys()):
-    print(val.evalf())
-# Eigen Values
-
## 3.14160000000000
-## 0.333362500589654
-## 2.71823749941035
-
pprint(eVects)
-# Eigen Vectors
-
##               [1.0]                            [0.666693567376375]            
-##               [   ]                            [                 ]            
-## [(3.1416, 1, [[1.0]]), (0.333362500589654, 1, [[3.66681936211085 ]]), (2.71823
-##               [   ]                            [                 ]            
-##               [ 0 ]                            [       1.0       ]            
-## 
-##                [-0.666648265089379]   
-##                [                  ]   
-## 749941035, 1, [[-5.66677405982385 ]])]
-##                [                  ]   
-##                [       1.0        ]
-
for val in list(eVals.keys()):
-    print("----------------------------------------------------------")
-    eig = val.evalf()
-    m = A-I*eig
-    print(eig)
-    pprint(m)
-    pprint(m.rref())
-# Eigen Value Matricies
-
## ----------------------------------------------------------
-## 3.14160000000000
-## [0.2693  -0.2693  -1.0643]
-## [                        ]
-## [1.587   -1.587   -5.5361]
-## [                        ]
-## [0.2981  -0.2981  -1.9139]
-##  [1  0  0]            
-##  [       ]            
-## ([0  1  0], (0, 1, 2))
-##  [       ]            
-##  [0  0  1]            
-## ----------------------------------------------------------
-## 0.333362500589654
-## [3.07753749941035      -0.2693            -1.0643     ]
-## [                                                     ]
-## [     1.587        1.22123749941035       -5.5361     ]
-## [                                                     ]
-## [     0.2981           -0.2981       0.894337499410346]
-##  [1  0  -0.666693567376375]         
-##  [                        ]         
-## ([0  1  -3.66681936211085 ], (0, 1))
-##  [                        ]         
-##  [0  0          0         ]         
-## ----------------------------------------------------------
-## 2.71823749941035
-## [0.692662500589653       -0.2693            -1.0643     ]
-## [                                                       ]
-## [      1.587        -1.16363749941035       -5.5361     ]
-## [                                                       ]
-## [     0.2981             -0.2981       -1.49053749941035]
-##  [1  0  0]            
-##  [       ]            
-## ([0  1  0], (0, 1, 2))
-##  [       ]            
-##  [0  0  1]
-
- - - - -
- - - - - - - - diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/HW/1/eigen.py b/MATH5620_NumericalSolutionsOfDifferentialEquations/HW/1/eigen.py deleted file mode 100644 index 40e9202..0000000 --- a/MATH5620_NumericalSolutionsOfDifferentialEquations/HW/1/eigen.py +++ /dev/null @@ -1,25 +0,0 @@ -#!/usr/bin/env python - -from sympy import * - -init_printing() - -A = Matrix([[3.4109, -0.2693, -1.0643],[1.5870, 1.5546, -5.5361],[0.2981, -0.2981, 1.2277]]) -I = Matrix([[1,0,0],[0,1,0],[0,0,1]]) -y = symbol('y') -lamba = symbol('lamba') - -pprint(I) - -eVals = A.eigenvals() - -eVects = A.eigenvects() - -#print(eVals) -pprint(A) -print("\nEigen Values:") -for val in list(eVals.keys()): - print(val.evalf()) - -print("\nEigen Vectors:") -pprint(eVects) diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/HW/1/eigen.pyc b/MATH5620_NumericalSolutionsOfDifferentialEquations/HW/1/eigen.pyc deleted file mode 100644 index c72c1b1..0000000 Binary files a/MATH5620_NumericalSolutionsOfDifferentialEquations/HW/1/eigen.pyc and /dev/null differ diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/HW/1/eigen.rmd b/MATH5620_NumericalSolutionsOfDifferentialEquations/HW/1/eigen.rmd deleted file mode 100644 index 7152904..0000000 --- a/MATH5620_NumericalSolutionsOfDifferentialEquations/HW/1/eigen.rmd +++ /dev/null @@ -1,36 +0,0 @@ ---- -title: HW1.10 -author: Philip ---- - -#### Compute the eigenvalues and eigenvectors for the following matrix using ye olde proverbally pencil and paper - -```{python} -from sympy import * - -A = Matrix([[3.4109, -0.2693, -1.0643],[1.5870, 1.5546, -5.5361],[0.2981, -0.2981, 1.2277]]) -I = Matrix([[1,0,0],[0,1,0],[0,0,1]]) - -eVals = A.eigenvals() -eVects = A.eigenvects() - -pprint(A) -# A matrix - -for val in list(eVals.keys()): - print(val.evalf()) -# Eigen Values - -pprint(eVects) -# Eigen Vectors - -for val in list(eVals.keys()): - print("----------------------------------------------------------") - eig = val.evalf() - m = A-I*eig - print(eig) - pprint(m) - pprint(m.rref()) -# Eigen Value Matricies - -``` diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/SoftwareManual.md b/MATH5620_NumericalSolutionsOfDifferentialEquations/SoftwareManual.md index 5a1438e..e17c157 100644 --- a/MATH5620_NumericalSolutionsOfDifferentialEquations/SoftwareManual.md +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/SoftwareManual.md @@ -1,11 +1,163 @@ --- +layout: default title: MATH 5620 Software Manual --- + Home + # Table of Contents +_(Homework problems below)_ + ### Basic Routines -- [Machine Epsilon](https://philipnelson5.github.io/class-projects/MATH5620_NumericalSolutionsOfDifferentialEquations/machineEpsilon/manual) -- [Absolute and Relative Error](https://philipnelson5.github.io/class-projects/MATH5620_NumericalSolutionsOfDifferentialEquations/error/manual) -- [Logistic Differential Equation](https://philipnelson5.github.io/class-projects/MATH5620_NumericalSolutionsOfDifferentialEquations/logistic/manual) -- [Second Order Linear DE with Constant Coefficents](https://philipnelson5.github.io/class-projects/MATH5620_NumericalSolutionsOfDifferentialEquations/secondOrderLinear/manual) +- [Machine Epsilon](./machineEpsilon/manual) +- [Absolute and Relative Error](./error/manual) +- [Logistic Differential Equation](./logistic/manual) +- [Second Order Linear DE with Constant Coefficients](./secondOrderLinear/manual) + +### Norms +- [Vector P Norm](./matrix/manual_vector_pnorm) +- [Vector Infinity Norm](./matrix/manual_vector_infinity_norm) +- [Matrix One Norm](./matrix/manual_matrix_one_norm) +- [Matrix Infinity Norm](./matrix/manual_matrix_infinity_norm) + +### Linear Solvers +- [Thomas Algorithm](./matrix/manual_thomas_algorithm) +- [Jacobi Iteration](./matrix/manual_jacobi_iteration) +- [Conjugate Gradient](./conjugateGradient/manual_conjugate_gradient) +- [Linear Solver by LU Factorization](./matrix/manual_linear_solve_lu) + - [Lu Factorization](./matrix./manual_lu_factorization) + - [Forward Substitution](./matrix./manual_forward_sub) + - [Back Substitution](./matrix./manual_back_sub) + +### Eigenvalue Methods +- [Power Method Iteration for Largest Eigenvalue](./matrix/manual_power_iteration) +- [Inverse Power Method for Smallest Eigenvalue](./matrix/manual_inverse_power_iteration) +- [Power Method Iteration for Solving 2nd Order FD Elliptic ODE](./matrix/example_power_iteration_elliptic_ode) +- [Inverse Power Method for Solving 2nd Order FD Elliptic ODE](./matrix/example_inverse_power_iteration_elliptic_ode) + +### Elliptic Problems +- [Finite Difference Coefficients](./finiteDiffMethods/manual_finite_diff_coeff) +- [Initialize Elliptic ODE](./finiteDiffMethods/manual_init_elliptic_ode) +- [Solve Elliptic ODE](./finiteDiffMethods/manual_solve_elliptic_ode) +- [Solve Laplace Equation with 5-point Stencil](./matrix/manual_solve_five_point_stencil) + - [Generate 5-point Stencil](./matrix/manual_gen_five_point_stencil) + - [Generate Mesh](./matrix/manual_gen_mesh) + - [Initialize B for Mesh](./matrix/manual_init_b) + +### First Order IVPs +- [Explicit Euler](./explicitEuler/manual_explicit_euler) +- [Implicit Euler](./implicitEuler/manual_implicit_euler) + - [Newton's Method](./newtonsMethod/manual_newtons_method) +- [Runge Kutta order 2](./rungeKuttaOrder2/manual_runge_kutta_order2) +- [Runge Kutta order 4](./rungeKuttaOrder4/manual_runge_kutta_order4) +- [Predicticor Corrector - Adams Bashforth / Adams Moulton](./predictorCorrector/manual_predictor_corrector) + +### Parabolic Problems +- [Upwinding](./upwinding/manual_upwinding) +- [Lax-Wendorff Method](./laxWendroff/manual_lax_wendroff) +- [Warming and Beam Method](./warmingAndBeam/manual_warming_and_beam) + +--- + +### Homework 1 +*due: 25 January 2018* + +| Problem | Software Manual| +| :-----------------|:---------------| +| **Problem 2-3.** | [Machine Epsilon](./machineEpsilon/manual)| +| **Problem 4.** | [Absolute and Relative Error](./error/manual)| +| **Problem 6.** | [Logistic Differential Equation](./logistic/manual)| +| **Problem 7.** | [Second Order Linear DE with Constant Coefficients](./secondOrderLinear/manual)| + + +### Homework 2 +*due: 8 February 2018* + +| Problem | Software Manual| +| :-----------------|:---------------| +| **Problem 1-2.** | [Finite Difference Coefficients](./finiteDiffMethods/manual_finite_diff_coeff)| +| **Problem 3.** | [Initialize Elliptic ODE](./finiteDiffMethods/manual_init_elliptic_ode)| +| **Problem 4.** | [Thomas Algorithm](./matrix/manual_thomas_algorithm)| +| **Problem 5.** | [Linear Solver by LU Factorization](./matrix/manual_linear_solve_lu)| +| | [Lu Factorization](./matrix/manual_lu_factorize)| +| | [Forward Substitution](./matrix/manual_forward_sub)| +| | [Back Substitution](./matrix/manual_back_sub)| +| **Problem 6.** | [Jacobi Iteration](./matrix/manual_jacobi_iteration)| +| **Problem 8.** | [Solve Elliptic ODE](./finiteDiffMethods/manual_solve_elliptic_ode)| +| | [Vector P Norm](./matrix/manual_vector_pnorm)| +| | [Vector Infinity Norm](./matrix/manual_vector_infinity_norm)| + +### Homework 3 +*due: 1 March 2018* + +| Problem | Software Manual| +| :-----------------|:---------------| +| **Problem 1.** | [Matrix One Norm](./matrix/manual_matrix_one_norm)| +| | [Matrix Infinity Norm](./matrix/manual_matrix_infinity_norm)| +| **Problem 2.** | [Power Method Iteration for Largest Eigenvalue](./matrix/manual_power_iteration)| +| | [Inverse Power Method for Smallest Eigenvalue](./matrix/manual_inverse_power_iteration)| +| **Problem 3.** | [Power Method Iteration for Solving 2nd Order FD Elliptic ODE](./matrix/example_power_iteration_elliptic_ode)| +| **Problem 4.** | [Inverse Power Method for Solving 2nd Order FD Elliptic ODE](./matrix/example_inverse_power_iteration_elliptic_ode)| +| **Problem 5.** | [Solve Laplace Equation with 5-point Stencil](./matrix/manual_solve_five_point_stencil)| +| | [Generate 5-point Stencil](./matrix/manual_gen_five_point_stencil)| +| | [Generate Mesh](./matrix/manual_gen_mesh)| +| | [Initialize B for Mesh](./matrix/manual_init_b)| +| **Problem 6.** | [Solve Laplace Equation with 9-point Stencil](./matrix/manual_solve_nine_point_stencil)| +| | [Generate 9-point Stencil](./matrix/manual_gen_nine_point_stencil)| + +### Homework 4 +*due: 22 March 2018* + +| Problem | Software Manual| +| :-----------------|:---------------| +| **Problem 1.** | [Gauss-Seidel](./gaussSidel/manual_gauss_sidel)| +| **Problem 2.** | [Conjugate Gradient](./conjugateGradient/manual_conjugate_gradient)| +| **Problem 3.** | [Application of Conjugate Gradient](./testConjugateGradientFivePoint/manual_solve_five_point_stencil_test)| +| **Problem 4.** | [Explicit Euler](./explicitEuler/manual_explicit_euler)| + +### Homework 5 +*due: 3 April 2018* + +| Problem | Software Manual| +| :-----------------|:---------------| +| **Problem 1.** | [First Order IVP Test](./5.1IVP/IVP_test)| +| | [Logistic Model Test](./logistic2/manual)| +| **Problem 2.** | [IVP via Explicit Euler Test](./explicitEulerTest/manual_explicit_euler_test)| +| **Problem 3.** | [Implicit Euler](./implicitEuler/manual_implicit_euler)| +| | [Newton's Method](./newtonsMethod/manual_newtons_method)| +| **Problem 4.** | [Runge Kutta order 2](./rungeKuttaOrder2/manual_runge_kutta_order2)| +| | [Runge Kutta order 4](./rungeKuttaOrder4/manual_runge_kutta_order4)| +| **Problem 5.** | [Adam's Bashford](./predictorCorrector/manual_predictor_corrector)| +| **Problem 6.** | [Summary of Iterative Methods](./SumaryOfIterativeMethods)| + +### Homework 6 +*due: 24 April 2018* + +| Problem | Software Manual| +| :-----------------|:---------------| +| **Write up** | [Experiments 7.1, 7.2, 7.4](./hw6_experiments.md)| + +### Homework 7 +*due: 5 May 2018* + +| Problem | Software Manual| +| :-----------------|:---------------| +| **Problem 1.** | [Heat Equation - Explicit Euler](./heatEquations/manual_heat_equation_explicit_euler)| +| **Problem 2.** | [Heat Equation - Implicit Euler](./heatEquations/manual_heat_equation_implicit_euler)| +| **Problem 4.** | [Heat Equation - Predictor Corrector](./heatEquations/manual_heat_equation_predictor_corrector)| +| **Problem 5.** | [Heat Equation - Runge Kutta Order 4](./heatEquations/manual_heat_equation_runge_kutta)| +{% comment %} +| **Problem 3.** | [Changing Time Time Step](./heatEquations/manual_heat_equation)| +{% endcomment %} + +### Homework 8 +*due: 5 May 2018* + +| Problem | Software Manual| +| :-----------------|:---------------| +| **Problem 1.1** | [Upwinding](./upwinding/manual_upwinding)| +| **Problem 1.2** | [Lax-Wendorff Method](./laxWendroff/manual_lax_wendroff)| +| **Problem 1.3** | [Warming and Beam Method](./warmingAndBeam/manual_warming_and_beam)| +| **Problem 2.1** | [vonNeuman Stability Analysis Lax Wendroff](./stabilityAnalysis/laxWendroff)| +| **Problem 2.2** | [vonNeuman Stability Analysis Warming and Beam](./stabilityAnalysis/warmingAndBeam)| diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/SumaryOfIterativeMethods.md b/MATH5620_NumericalSolutionsOfDifferentialEquations/SumaryOfIterativeMethods.md new file mode 100644 index 0000000..e9cf6c8 --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/SumaryOfIterativeMethods.md @@ -0,0 +1,19 @@ +--- +title: IVP Test +math: true +layout: default +--- + +{% include mathjax.html %} + +# Iterative Methods Summary + +Iterative methods studied in assignment 5 are examples iterative methods used to solve initial value problems. +Each method has it's strengths and weaknesses, there is best method. Each has different conditions for convergence, but some have better overall performance. + +The Explicit Euler Method is a very simple technique compared to Runge Kutta and Adams Bashforth and Adams Moulton. This method represents a single step of integration at a point. The Implicit Euler technique, however uses the unknown next value of the iteration in the same setup. This results in an equation that must be solved using some other method such as the Newton method for finding zeroes. The Implicit Euler Method has order one. This means that the local truncation error is \\(O(h^{2})\\). The error at a specific time \\(t\\) is \\(O(h)\\). + +Runge Kutta, which is favored of engineers, is favored due to its stability and ease of use. Runge Kutta techniques exist for arbitrary order, the Kunge Kutta order 4 is the most widely known and is referred to as the "classical Runge Kutta method". + +The Predictor Corrector method used Adams Bashforth for the prediction step then Adams Moulton for the corrector. This is the typical form that predictor-corrector methods take, an explicit method for the predictor step and an implicit method for the corrector step. The Adams–Bashforth methods are explicit methods. The coefficients are \\(a_{s-1}=-1\\) and \\(a_{s-2}=\cdots =a_0=0\\), while the \\(b_j\\) are chosen such that the methods has order s (this determines the methods uniquely). The Adams–Moulton methods are similar to the Adams–Bashforth methods in that they also have \\(a_{s-1}=-1\\) and \\(a_{s-2}=\cdots =a_0=0\\). Again the b coefficients are chosen to obtain the highest order possible. However, the Adams–Moulton methods are implicit methods. By removing the restriction that \\(b_s = 0\\) , an s-step Adams–Moulton method can reach order s+1, while an s-step Adams–Bashforth methods has only order s. + diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/_config.yml b/MATH5620_NumericalSolutionsOfDifferentialEquations/_config.yml new file mode 100644 index 0000000..060290a --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/_config.yml @@ -0,0 +1,3 @@ +theme: jekyll-theme-modernist +markdown: kramdown +mathjax: true diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/_includes/mathjax.html b/MATH5620_NumericalSolutionsOfDifferentialEquations/_includes/mathjax.html new file mode 100644 index 0000000..f7c4557 --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/_includes/mathjax.html @@ -0,0 +1,3 @@ + diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/_layouts/default.html b/MATH5620_NumericalSolutionsOfDifferentialEquations/_layouts/default.html new file mode 100644 index 0000000..3866843 --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/_layouts/default.html @@ -0,0 +1,63 @@ + + + + + + +{% seo %} + + + + + + + + + +
+
+

{{ site.title | default: site.github.repository_name }}

+ {% if site.description or site.github.project_tagline %} +

{{ site.description | default: site.github.project_tagline }}

+ {% endif %} + + +
+
+ + {{ content }} + +
+
+
+ {% if site.github.is_project_page %} + + {% endif %} + +
+ + + {% if site.google_analytics %} + + {% endif %} + + diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/assets/css/style.scss b/MATH5620_NumericalSolutionsOfDifferentialEquations/assets/css/style.scss new file mode 100644 index 0000000..8f64127 --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/assets/css/style.scss @@ -0,0 +1,8 @@ +--- +--- + +@import "{{ site.theme }}"; + +.wrapper { + width: 800px; +} diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/conjugateGradient/Makefile b/MATH5620_NumericalSolutionsOfDifferentialEquations/conjugateGradient/Makefile new file mode 100644 index 0000000..01ada52 --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/conjugateGradient/Makefile @@ -0,0 +1,13 @@ +OBJS = ../matrix/matrix.hpp ../matrix/matrix_util.hpp ../matrix/vector_util.hpp ../gaussSidel/gaussSidel.hpp conjugateGradient.hpp main.cpp + +CC = g++ +FLAGS = -std=c++17 +DEBUG_FLAGS = -g3 -Og +RELEASE_FLAGS = -g0 -O3 +WARNINGS = -Wall -Wextra -Werror + +release: $(OBJS) + $(CC) $(FLAGS) $(RELEASE_FLAGS) $(WARNINGS) $(OBJS) -o conjugateGradient.out + +debug: $(OBJS) + $(CC) $(FLAGS) $(DEBUG_FLAGS) $(WARNINGS) $(OBJS) -o debug.out diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/conjugateGradient/conjugateGradient.hpp b/MATH5620_NumericalSolutionsOfDifferentialEquations/conjugateGradient/conjugateGradient.hpp new file mode 100644 index 0000000..95166f7 --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/conjugateGradient/conjugateGradient.hpp @@ -0,0 +1,49 @@ +#ifndef CONJUGATE_GRADIENT_HPP +#define CONJUGATE_GRADIENT_HPP + +#include "../machineEpsilon/maceps.hpp" +#include "../matrix/matrix.hpp" +#include "../matrix/vector_util.hpp" +#include + +template +std::array conjugate_gradient(Matrix& A, + std::array const& b, + unsigned int const& MAX_ITERATIONS = 1000u) +{ + auto ct = 0u; + auto tol = maceps().maceps; + std::array x_k, x_k1; + x_k.fill(0); + x_k1.fill(0); + auto r_k = b; + auto r_k1 = r_k, r_k2 = r_k; + auto p_k = r_k, p_k1 = r_k; + for (auto k = 0u; k < MAX_ITERATIONS; ++k) + { + ++ct; + if (k != 0) + { + auto b_k = (r_k1 * r_k1) / (r_k2 * r_k2); + p_k = r_k1 + b_k * p_k1; + } + auto s_k = A * p_k; + auto a_k = r_k1 * r_k1 / (p_k * s_k); + x_k = x_k1 + a_k * p_k; + r_k = r_k1 - a_k * s_k; + + if (allclose(x_k, x_k1, tol)) + { + std::cout << "Conjugate Gradient completed in " << ct << " iterations\n"; + return x_k; + } + + r_k2 = r_k1; + r_k1 = r_k; + x_k1 = x_k; + p_k1 = p_k; + } + return x_k; +} + +#endif diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/conjugateGradient/main.cpp b/MATH5620_NumericalSolutionsOfDifferentialEquations/conjugateGradient/main.cpp new file mode 100644 index 0000000..3b91fb3 --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/conjugateGradient/main.cpp @@ -0,0 +1,21 @@ +#include "../gaussSidel/gaussSidel.hpp" +#include "../matrix/matrix.hpp" +#include "conjugateGradient.hpp" +#include + +int main() +{ + Matrix A({{16, 3}, {7, -11}}); + std::array b({{11, 13}}); + std::cout << "A\n" << A << "b\n" << b << '\n'; + std::cout << gauss_sidel(A, b, 100u) << "\n"; + std::cout << A.jacobiIteration(b, 100u) << "\n"; + std::cout << conjugate_gradient(A, b, 1000u) << "\n"; + + Matrix A1({{10, -1, 2, 0}, {-1, 11, -1, 3}, {2, -1, 10, -1}, {0, 3, -1, 8}}); + std::array b1({{6, 25, -11, 15}}); + std::cout << "A\n" << A1 << "b\n" << b1 << '\n'; + std::cout << gauss_sidel(A1, b1, 100u) << "\n"; + std::cout << A1.jacobiIteration(b1, 100u) << "\n"; + std::cout << conjugate_gradient(A1, b1, 1000u) << "\n"; +} diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/conjugateGradient/manual_conjugate_gradient.md b/MATH5620_NumericalSolutionsOfDifferentialEquations/conjugateGradient/manual_conjugate_gradient.md new file mode 100644 index 0000000..4a8f58a --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/conjugateGradient/manual_conjugate_gradient.md @@ -0,0 +1,151 @@ +--- +title: Conjugate Gradient +math: true +layout: default +--- + +{% include mathjax.html %} + + Table of Contents +# Conjugate Gradient Method + +**Routine Name:** `conjugate_gradient` + +**Author:** Philip Nelson + +**Language:** C++ + +## Description + +`conjugate_gradient` is an algorithm for the numerical solution of particular systems of linear equations, namely those whose matrix is symmetric and positive-definite. `conjugate_gradient` is implemented as an iterative algorithm, applicable to sparse systems that are too large to be handled by a direct implementation. [1](https://en.wikipedia.org/wiki/Conjugate_gradient_method) + +## Input + + +``` +std::array conjugate_gradient(Matrix& A, + std::array const& b, + unsigned int const& MAX_ITERATIONS = 1000u) +``` +requires: + +* `Matrix& A` - an `N`x`N` matrix of type `T` +* `std::array const& b` - the b vector +* `unsigned int const& MAX_ITERATIONS` - the maximum iterations (default 1000) + +## Output + +The solution vector \\(x\\). + +## Code +{% highlight c++ %} +template +std::array conjugate_gradient(Matrix& A, + std::array const& b, + unsigned int const& MAX_ITERATIONS = 1000u) +{ + auto ct = 0u; + auto tol = maceps().maceps; + std::array x_k, x_k1; + x_k.fill(0); + x_k1.fill(0); + auto r_k = b; + auto r_k1 = r_k, r_k2 = r_k; + auto p_k = r_k, p_k1 = r_k; + for (auto k = 0u; k < MAX_ITERATIONS; ++k) + { + ++ct; + if (k != 0) + { + auto b_k = (r_k1 * r_k1) / (r_k2 * r_k2); + p_k = r_k1 + b_k * p_k1; + } + auto s_k = A * p_k; + auto a_k = r_k1 * r_k1 / (p_k * s_k); + x_k = x_k1 + a_k * p_k; + r_k = r_k1 - a_k * s_k; + + if (allclose(x_k, x_k1, tol)) + { + std::cout << "Conjugate Gradient completed in " << ct << " iterations\n"; + return x_k; + } + + r_k2 = r_k1; + r_k1 = r_k; + x_k1 = x_k; + p_k1 = p_k; + } + return x_k; +} +{% endhighlight %} + +## Example +{% highlight c++ %} +int main() +{ + Matrix A({ + {16, 3}, + { 7, -11} + }); + std::array b({ + {11, 13} + }); + std::cout << "A\n" << A << "b\n" << b << '\n'; + std::cout << gauss_sidel(A, b, 100u) << "\n"; + std::cout << A.jacobiIteration(b, 100u) << "\n"; + std::cout << conjugate_gradient(A, b, 10000u) << "\n"; + + Matrix A1({ + {10, -1, 2, 0}, + {-1, 11, -1, 3}, + { 2, -1, 10, -1}, + { 0, 3, -1, 8} + }); + std::array b1({ + {6, 25, -11, 15} + }); + std::cout << "A\n" << A1 << "b\n" << b1 << '\n'; + std::cout << gauss_sidel(A1, b1, 100u) << "\n"; + std::cout << A1.jacobiIteration(b1, 100u) << "\n"; + std::cout << conjugate_gradient(A1, b1, 1000u) << "\n"; +} +{% endhighlight %} + +## Result +As opposed to Gauss Sidel, Conjugate Gradient converges differently on the solution depending on the topology of the region around the solution. +``` +A +| 16 3 | +| 7 -11 | +b +[ 11 13 ] + +Gauss Sidel completed in 18 iterations +[ 0.812 -0.665 ] + +Jacobi Iteration completed in 35 iterations +[ 0.812 -0.665 ] + +Conjugate Gradient completed in 75 iterations +[ 0.812 -0.665 ] + +A +| 10 -1 2 0 | +| -1 11 -1 3 | +| 2 -1 10 -1 | +| 0 3 -1 8 | +b +[ 6 25 -11 15 ] + +Gauss Sidel completed in 17 iterations +[ 1 2 -1 1 ] + +Jacobi Iteration completed in 43 iterations +[ 1 2 -1 1 ] + +Conjugate Gradient completed in 5 iterations +[ 1 2 -1 1 ] +``` + +**Last Modification date:** 21 March 2018 diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/error/main.cpp b/MATH5620_NumericalSolutionsOfDifferentialEquations/error/main.cpp index 1c4ee2b..0f34fdb 100644 --- a/MATH5620_NumericalSolutionsOfDifferentialEquations/error/main.cpp +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/error/main.cpp @@ -1,13 +1,14 @@ #include "error.hpp" -#include #include +#include int main() { - double approx = 3.2; - double value = 3.14159; + auto approx = 3.1416; + auto value = M_PI; - std::cout << "Approximate: " << approx << "\nReal Value: " << value << std::endl + std::cout << "Approximate: " << approx << "\nReal Value: " << value + << std::endl << std::endl; std::cout << "Absolute: " << absoluteError(approx, value) << std::endl; std::cout << "Relative: " << relativeError(approx, value) << std::endl; diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/error/manual.md b/MATH5620_NumericalSolutionsOfDifferentialEquations/error/manual.md index c536028..3302e32 100644 --- a/MATH5620_NumericalSolutionsOfDifferentialEquations/error/manual.md +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/error/manual.md @@ -6,7 +6,7 @@ layout: default {% include mathjax.html %} - Table of Contents + Table of Contents # Error **Routine Names:** [absoluteError](#input-absoluteerror) and [relativeError](#input-relativeerror) diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/explicitEuler/Makefile b/MATH5620_NumericalSolutionsOfDifferentialEquations/explicitEuler/Makefile new file mode 100644 index 0000000..02cc602 --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/explicitEuler/Makefile @@ -0,0 +1,13 @@ +OBJS = ../matrix/matrix.hpp ../matrix/matrix_util.hpp ../matrix/vector_util.hpp main.cpp + +CC = g++ +FLAGS = -std=c++17 +DEBUG_FLAGS = -g3 -Og +RELEASE_FLAGS = -g0 -O3 +WARNINGS = -Wall -Wextra -Werror + +release: $(OBJS) + $(CC) $(FLAGS) $(RELEASE_FLAGS) $(WARNINGS) $(OBJS) -o explicitEuler.out + +debug: $(OBJS) + $(CC) $(FLAGS) $(DEBUG_FLAGS) $(WARNINGS) $(OBJS) -o debug.out diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/explicitEuler/explicitEuler.hpp b/MATH5620_NumericalSolutionsOfDifferentialEquations/explicitEuler/explicitEuler.hpp new file mode 100644 index 0000000..bd336d4 --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/explicitEuler/explicitEuler.hpp @@ -0,0 +1,20 @@ +#ifndef EXPLICIT_EULER_HPP +#define EXPLICIT_EULER_HPP + +#include "../machineEpsilon/maceps.hpp" +#include +#include + +template +T explicit_euler(T x0, T y0, T x, T dt, F f) +{ + auto tol = maceps().maceps; + while (std::abs(x - x0) > tol) + { + y0 = y0 + (dt * f(x0, y0)); + x0 += dt; + } + return y0; +} + +#endif diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/explicitEuler/main.cpp b/MATH5620_NumericalSolutionsOfDifferentialEquations/explicitEuler/main.cpp new file mode 100644 index 0000000..e97758e --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/explicitEuler/main.cpp @@ -0,0 +1,6 @@ +#include "explicitEuler.hpp" +#include + +int main() { + std::cout << explicit_euler(0.0, -1.0, 0.4, 0.1, [](double a, double b){return a*a+2*b;}) << std::endl; +} diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/explicitEuler/manual_explicit_euler.md b/MATH5620_NumericalSolutionsOfDifferentialEquations/explicitEuler/manual_explicit_euler.md new file mode 100644 index 0000000..4a6c90b --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/explicitEuler/manual_explicit_euler.md @@ -0,0 +1,67 @@ +--- +title: Explicit Euler +math: true +layout: default +--- + +{% include mathjax.html %} + + Table of Contents +# Explicit Euler Method + +**Routine Name:** `explicit_euler` + +**Author:** Philip Nelson + +**Language:** C++ + +## Description + +`explicit_euler` is a first-order numerical procedure for solving ordinary differential equations (ODEs) with a given initial value. It is the most basic explicit method for numerical integration of ordinary differential equations and is the simplest Runge–Kutta method. [1](https://en.wikipedia.org/wiki/Euler_method) + +## Input + + +`explicit_euler(T x0, T y0, T x, T dt, F f)` requires: + +* `T x0` - the initial `x` +* `T y0` - the initial `y` +* `T x` - the value of `x` for which you want to find value of `y` +* `T dt` - the delta t step +* `F f` - the function defining \\(\frac{dy}{dx}\\) + +## Output + +The value of `y` at `x`. + +## Code +{% highlight c++ %} +template +T explicit_euler(T x0, T y0, T x, T dt, F f) +{ + auto tol = maceps().maceps; + while (std::abs(x - x0) > tol) + { + y0 = y0 + (dt * f(x0, y0)); + x0 += dt; + } + return y0; +} +{% endhighlight %} + +## Example +{% highlight c++ %} +int main() +{ + std::cout << + explicit_euler(0.0, -1.0, 0.4, 0.1, [](double a, double b){return a*a+2*b;}) + << std::endl; +} +{% endhighlight %} + +## Result +``` +-2.05836 +``` + +**Last Modification date:** 21 March 2018 diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/explicitEulerTest/Makefile b/MATH5620_NumericalSolutionsOfDifferentialEquations/explicitEulerTest/Makefile new file mode 100644 index 0000000..1057e95 --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/explicitEulerTest/Makefile @@ -0,0 +1,13 @@ +OBJS = ../matrix/matrix.hpp ../matrix/matrix_util.hpp ../matrix/vector_util.hpp main.cpp + +CC = g++ +FLAGS = -std=c++17 +DEBUG_FLAGS = -g3 -Og +RELEASE_FLAGS = -g0 -O3 +WARNINGS = -Wall -Wextra -Werror + +release: $(OBJS) + $(CC) $(FLAGS) $(RELEASE_FLAGS) $(WARNINGS) $(OBJS) -o explicitEulerTest.out + +debug: $(OBJS) + $(CC) $(FLAGS) $(DEBUG_FLAGS) $(WARNINGS) $(OBJS) -o debug.out diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/explicitEulerTest/main.cpp b/MATH5620_NumericalSolutionsOfDifferentialEquations/explicitEulerTest/main.cpp new file mode 100644 index 0000000..57af737 --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/explicitEulerTest/main.cpp @@ -0,0 +1,45 @@ +#include "../logistic/logistic.hpp" +#include "../explicitEuler/explicitEuler.hpp" +#include "../5.1IVP/firstOrderIVP.hpp" +#include +#include + +int main() { + double dt = 0.00001; + + double alpha = 10.0; + std::vector lambdas = { 1.0, -1.0, 100.0 }; + + double gamma = 0.1; + double beta = 0.0001; + std::vector Pos = { 25.0, 40000.0 }; + + std::cout << "----- Lambda Differential Equation -----" << std::endl; + for (const auto lambda : lambdas) { + std::cout << explicit_euler( alpha, beta, gamma, dt, + [=](double a, double b) { + (void)b; + return lambda * a; + } + ) << '\n'; + + auto solveIVP = firstOrderIVPSolver(lambda, alpha); + + std::cout << solveIVP(dt); + } + + std::cout << std::endl; + std::cout << "----- Logistic Differential Equation -----" << std::endl; + for (const auto p0 : Pos) { + std::cout << explicit_euler( alpha, beta, gamma, dt, + [=](double a, double b) { + (void)b; + return gamma * a - beta * a * a; + } + ) << '\n'; + + std::cout << logistic(alpha, beta, dt, p0) << '\n'; + } + + return EXIT_SUCCESS; +} diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/explicitEulerTest/manual_explicit_euler_test.md b/MATH5620_NumericalSolutionsOfDifferentialEquations/explicitEulerTest/manual_explicit_euler_test.md new file mode 100644 index 0000000..db370fc --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/explicitEulerTest/manual_explicit_euler_test.md @@ -0,0 +1,139 @@ +--- +title: Explicit Euler Tests +math: true +layout: default +--- + +{% include mathjax.html %} + + Table of Contents +# Explicit Euler Method Tests + +**Author:** Philip Nelson + +**Language:** C++ + +## Description + +This tests the Explicit Euler method against the test cases in problem 1 of this assignment + +## Code +{% highlight c++ %} +int main() +{ + double dt = 0.00001; + + double alpha = 10.0; + std::vector lambdas = { 1.0, -1.0, 100.0 }; + + double gamma = 0.1; + double beta = 0.0001; + std::vector Pos = { 25.0, 40000.0 }; + + std::cout << "----- Lambda Differential Equation -----" << std::endl; + for (const auto lambda : lambdas) { + std::cout << explicit_euler( alpha, beta, gamma, dt, + [=](double a, double b) { + (void)b; + return lambda * a; + } + ) << '\n'; + + auto solveIVP = firstOrderIVPSolver(lambda, alpha); + + std::cout << solveIVP(dt); + } + + std::cout << std::endl; + std::cout << "----- Logistic Differential Equation -----" << std::endl; + for (const auto p0 : Pos) { + std::cout << explicit_euler( alpha, beta, gamma, dt, + [=](double a, double b) { + (void)b; + return gamma * a - beta * a * a; + } + ) << '\n'; + + std::cout << logistic(alpha, beta, dt, p0) << '\n'; + } + + return EXIT_SUCCESS; +} +{% endhighlight %} + +## Result +``` +----- Lambda Differential Equation ----- +exact t = 0 -> 10 +approx t = 0 -> 10 +exact t = 0.2 -> 12.214 +approx t = 0.2 -> 12.214 +exact t = 0.4 -> 14.9182 +approx t = 0.4 -> 14.9182 +exact t = 0.6 -> 18.2212 +approx t = 0.6 -> 18.2211 +exact t = 0.8 -> 22.2554 +approx t = 0.8 -> 22.2553 +exact t = 1 -> 27.1828 +approx t = 1 -> 27.1824 + +lambda = -1 +exact t = 0 -> 10 +approx t = 0 -> 10 +exact t = 0.2 -> 8.18731 +approx t = 0.2 -> 8.1873 +exact t = 0.4 -> 6.7032 +approx t = 0.4 -> 6.70319 +exact t = 0.6 -> 5.48812 +approx t = 0.6 -> 5.4881 +exact t = 0.8 -> 4.49329 +approx t = 0.8 -> 4.49327 +exact t = 1 -> 3.67879 +approx t = 1 -> 3.67881 + +lambda = 100 +exact t = 0 -> 10 +approx t = 0 -> 10 +exact t = 0.2 -> 4.85165e+09 +approx t = 0.2 -> 4.80341e+09 +exact t = 0.4 -> 2.35385e+18 +approx t = 0.4 -> 2.30727e+18 +exact t = 0.6 -> 1.14201e+27 +approx t = 0.6 -> 1.10828e+27 +exact t = 0.8 -> 5.54062e+35 +approx t = 0.8 -> 5.32351e+35 +exact t = 1 -> 2.68812e+44 +approx t = 1 -> 2.55455e+44 + +----- Logistic Differential Equation ----- + +p0 = 25 +exact t = 0 -> 25 +approx t = 0 -> 25 +exact t = 0.2 -> 25.4922 +approx t = 0.2 -> 25.4922 +exact t = 0.4 -> 25.9937 +approx t = 0.4 -> 25.9937 +exact t = 0.6 -> 26.5049 +approx t = 0.6 -> 26.5049 +exact t = 0.8 -> 27.0259 +approx t = 0.8 -> 27.0259 +exact t = 1 -> 27.5568 +approx t = 1 -> 27.5568 + +p0 = 40000 +exact t = 0 -> 40000 +approx t = 0 -> 40000 +exact t = 0.2 -> 22570.2 +approx t = 0.2 -> 22569.9 +exact t = 0.4 -> 15815.2 +approx t = 0.4 -> 15815 +exact t = 0.6 -> 12228 +approx t = 0.6 -> 12227.8 +exact t = 0.8 -> 10003.8 +approx t = 0.8 -> 10003.7 +exact t = 1 -> 8490.15 +approx t = 1 -> 8490.11 +``` + +**Last Modification date:** 3 April 2018 diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/finDiffCoeff/Makefile b/MATH5620_NumericalSolutionsOfDifferentialEquations/finDiffCoeff/Makefile new file mode 100644 index 0000000..dbb1257 --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/finDiffCoeff/Makefile @@ -0,0 +1,12 @@ +OBJS = main.cpp finDiffCoeff.hpp + +CC = g++ +FLAGS = -std=c++17 +DEBUG_FLAGS = -Og -g3 -Wall -Wextra -Werror +RELEASE_FLAGS = -O3 + +release: $(OBJS) + $(CC) $(RELEASE_FLAGS) $(FLAGS) $(OBJS) -o finDiffCoeff.out + +debug: $(OBJS) + $(CC) $(DEBUG_FLAGS) $(FLAGS) $(OBJS) -o debug.out diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/finDiffCoeff/finDiffCoeff.hpp b/MATH5620_NumericalSolutionsOfDifferentialEquations/finDiffCoeff/finDiffCoeff.hpp new file mode 100644 index 0000000..8f4f6fd --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/finDiffCoeff/finDiffCoeff.hpp @@ -0,0 +1,48 @@ +#ifndef FIN_DIFF_COEFF_HPP +#define FIN_DIFF_COEFF_HPP + +#include "../matrix/matrix.hpp" +#include +#include + +unsigned int fact(unsigned int const& n) +{ + if (n == 0) + return 1; + unsigned int f = 1u; + for (auto i = 1u; i <= n; ++i) + f *= i; + return f; +} + +unsigned int binCoeff(unsigned int const& n, unsigned int const& k) +{ + auto nf = fact(n); + auto kf = fact(k); + auto nkf = fact(n - k); + + return nf / (kf * nkf); +} + +template +auto centralFinDiffCoeff() +{ + constexpr int size = 2.0 * std::floor((ord + 1.0) / 2.0) - 1.0 + acc; + constexpr int P = (size - 1.0) / 2.0; + + Matrix mat; + for (auto i = 0; i < size; ++i) + { + for (auto j = 0; j < size ; ++j) + { + mat[i][j] = std::pow(-P+j, i); + } + } + + std::array b; + b.fill(0.0); + b[ord] = fact(ord); + + return mat.solveLinearSystemLU(b); +} +#endif diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/finDiffCoeff/main.cpp b/MATH5620_NumericalSolutionsOfDifferentialEquations/finDiffCoeff/main.cpp new file mode 100644 index 0000000..6e6c8f2 --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/finDiffCoeff/main.cpp @@ -0,0 +1,15 @@ +#include "finDiffCoeff.hpp" +#include +#include "../matrix/matrix.hpp" +#include "../matrix/matrix_util.hpp" +#include "../matrix/vector_util.hpp" + +int main() +{ + auto coeffs = centralFinDiffCoeff(); + + std::cout << "coefficients of a second order derivative with 4th accuracy\n"; + std::cout << coeffs << std::endl; + + return EXIT_SUCCESS; +} diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/finiteDiffMethods/Makefile b/MATH5620_NumericalSolutionsOfDifferentialEquations/finiteDiffMethods/Makefile new file mode 100644 index 0000000..dbb1257 --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/finiteDiffMethods/Makefile @@ -0,0 +1,12 @@ +OBJS = main.cpp finDiffCoeff.hpp + +CC = g++ +FLAGS = -std=c++17 +DEBUG_FLAGS = -Og -g3 -Wall -Wextra -Werror +RELEASE_FLAGS = -O3 + +release: $(OBJS) + $(CC) $(RELEASE_FLAGS) $(FLAGS) $(OBJS) -o finDiffCoeff.out + +debug: $(OBJS) + $(CC) $(DEBUG_FLAGS) $(FLAGS) $(OBJS) -o debug.out diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/finiteDiffMethods/ellipticODE.hpp b/MATH5620_NumericalSolutionsOfDifferentialEquations/finiteDiffMethods/ellipticODE.hpp new file mode 100644 index 0000000..e03577b --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/finiteDiffMethods/ellipticODE.hpp @@ -0,0 +1,55 @@ +#ifndef ELLIPTIC_ODE_HPP +#define ELLIPTIC_ODE_HPP + +#include "finDiffCoeff.hpp" +#include + +template +std:: + tuple, std::array, std::array, std::array> + initEllipticODE(F f, T a, T b, T ua, T ub) +{ + auto h = (b - a) / N; + auto h2 = h * h; + + auto coeff = centralFinDiffCoeff(); + + std::array av, bv, cv, fv; + + av.fill(coeff[0]); + bv.fill(coeff[1]); + cv.fill(coeff[2]); + + for (auto i = 1u; i < N; ++i) + { + fv[i - 1] = h2 * f(a + i * h); + } + + fv[0] -= ua; + fv[fv.size() - 1] -= ub; + + return {av, bv, cv, fv}; +} + +template +std::array solveEllipticODE(F f, T a, T b, T ua, T ub) +{ + auto[_a, _b, _c, _f] = initEllipticODE(f, a, b, ua, ub); + + return Matrix::triDiagThomas(_a, _b, _c, _f); +} + + // template + // std::array solveEllipticK(F f, T a, T b, T ua, T ub, std::array k) + // { + // auto h = (b - a) / N; + // auto h2 = h * h; + // + // for (auto i = 1u; i < N; ++i) + // { + // fv[i - 1] = h2 * f(a + i * h); + // } + // + // } + +#endif diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/finiteDiffMethods/finDiffCoeff.hpp b/MATH5620_NumericalSolutionsOfDifferentialEquations/finiteDiffMethods/finDiffCoeff.hpp new file mode 100644 index 0000000..b2784a6 --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/finiteDiffMethods/finDiffCoeff.hpp @@ -0,0 +1,60 @@ +#ifndef FIN_DIFF_COEFF_HPP +#define FIN_DIFF_COEFF_HPP + +#include "../matrix/matrix.hpp" +#include +#include + +unsigned int fact(unsigned int const& n) +{ + if (n == 0) + return 1; + unsigned int f = 1u; + for (auto i = 1u; i <= n; ++i) + f *= i; + return f; +} + +unsigned int binCoeff(unsigned int const& n, unsigned int const& k) +{ + auto nf = fact(n); + auto kf = fact(k); + auto nkf = fact(n - k); + + return nf / (kf * nkf); +} + +template +auto centralFinDiffCoeff() +{ + constexpr int size = 2.0 * std::floor((ord + 1.0) / 2.0) - 1.0 + acc; + constexpr int P = (size - 1.0) / 2.0; + + Matrix mat; + for (auto i = 0; i < size; ++i) + { + for (auto j = 0; j < size; ++j) + { + mat[i][j] = std::pow(-P + j, i); + } + } + + std::array b; + b.fill(0.0); + b[ord] = fact(ord); + + return mat.solveLinearSystemLU(b); +} + +template +Matrix SecOrdFinDifMethEllipticMat() +{ + auto coeff = centralFinDiffCoeff(); + return Matrix ([&](int i, int j){ + if (i == j) return coeff[1]; + if (i+1 == j) return coeff[0]; + if (i == j+1) return coeff[2]; + return 0.0; + }); +} +#endif diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/finiteDiffMethods/main.cpp b/MATH5620_NumericalSolutionsOfDifferentialEquations/finiteDiffMethods/main.cpp new file mode 100644 index 0000000..821c65e --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/finiteDiffMethods/main.cpp @@ -0,0 +1,46 @@ +#include "../matrix/matrix.hpp" +#include "../matrix/matrix_util.hpp" +#include "../matrix/vector_util.hpp" +#include "ellipticODE.hpp" +#include "finDiffCoeff.hpp" +#include +#include + +double exact(double x) +{ + return 2.5 + 2.5 * x - 1 / (M_PI * M_PI) * std::sin(M_PI * x); +} + +int main() +{ + constexpr int N = 10; + auto a = 0.0, b = 1.0, ua = 2.5, ub = 5.0; + auto h = (b - a) / N; + + // auto[_a, _b, _c, _f] = initEllipticODE([](double x) { return std::sin(M_PI * x); }, a, + // b, ua, ub); + // + // std::cout << "a\n" << _a << std::endl; + // std::cout << "b\n" << _b << std::endl; + // std::cout << "c\n" << _c << std::endl; + // std::cout << "f\n" << _f << std::endl; + + auto calcVal = + solveEllipticODE([](double x) { return std::sin(M_PI * x); }, a, b, ua, ub); + + std::cout << calcVal << std::endl; + + std::array realVal; + for (auto i = 1u; i < N; ++i) + { + realVal[i - 1] = exact(a + i * h); + } + + std::cout << realVal << std::endl; + + std::cout << "P1 norm: " << pNorm(realVal - calcVal, 1) << std::endl; + std::cout << "P2 norm: " << pNorm(realVal - calcVal, 2) << std::endl; + std::cout << "inf norm: " << infNorm(realVal - calcVal) << std::endl; + + return EXIT_SUCCESS; +} diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/finiteDiffMethods/manual_elliptic_ode.md b/MATH5620_NumericalSolutionsOfDifferentialEquations/finiteDiffMethods/manual_elliptic_ode.md new file mode 100644 index 0000000..ba6e74c --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/finiteDiffMethods/manual_elliptic_ode.md @@ -0,0 +1,87 @@ +--- +title: Finite Difference Coefficients +math: true +layout: default +--- + +{% include mathjax.html %} + + Table of Contents +# Finite Difference Coefficients + +**Routine Name:** centralFinDiffCoeff + +**Author:** Philip Nelson + +**Language:** C++ + +## Description + +`centralFinDiffCoeff` returns a vector of the coefficients for finite difference approximations of an arbitrary order of accuracy for a given derivative. This routine uses binomial coefficients to calculate the coefficients with the formula below: + +``` +$ make +$ ./finDiffCoeff.out +``` + +This will compile and run the driver program. + +## Input + +`centralFinDiffCoeff()` takes no parameters, however it requires the type, order and accuracy as template parameters. + +* `T` - the type you want the coefficients in +* `ord` - the order of the derivative +* `acc` - the accuracy of the approximation + +## Output + +`centralFinDiffCoeff` returns a vector of coefficients. + +## Code +{% highlight c++ %} +template +auto centralFinDiffCoeff() +{ + constexpr int size = 2.0 * std::floor((ord + 1.0) / 2.0) - 1.0 + acc; + constexpr int P = (size - 1.0) / 2.0; + + Matrix mat; + for (auto i = 0; i < size; ++i) + { + for (auto j = 0; j < size ; ++j) + { + mat[i][j] = std::pow(-P+j, i); + } + } + + std::array b; + b.fill(0.0); + b[ord] = fact(ord); + + return mat.solveLinearSystemLU(b); +} +{% endhighlight %} + +## Example +{% highlight c++ %} +int main() +{ + auto coeffs = centralFinDiffCoeff(); + + std::cout << "coefficients of a second order derivative with 4th accuracy\n\n"; + std::cout << coeffs << std::endl; + + return EXIT_SUCCESS; +} +{% endhighlight %} + +## Result +``` +coefficients of a second order derivative with 4th accuracy + +[ -0.0833 1.33 -2.5 1.33 -0.0833 ] + +``` + +**Last Modification date:** 11 January 2018 diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/finiteDiffMethods/manual_finite_diff_coeff.md b/MATH5620_NumericalSolutionsOfDifferentialEquations/finiteDiffMethods/manual_finite_diff_coeff.md new file mode 100644 index 0000000..ba6e74c --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/finiteDiffMethods/manual_finite_diff_coeff.md @@ -0,0 +1,87 @@ +--- +title: Finite Difference Coefficients +math: true +layout: default +--- + +{% include mathjax.html %} + + Table of Contents +# Finite Difference Coefficients + +**Routine Name:** centralFinDiffCoeff + +**Author:** Philip Nelson + +**Language:** C++ + +## Description + +`centralFinDiffCoeff` returns a vector of the coefficients for finite difference approximations of an arbitrary order of accuracy for a given derivative. This routine uses binomial coefficients to calculate the coefficients with the formula below: + +``` +$ make +$ ./finDiffCoeff.out +``` + +This will compile and run the driver program. + +## Input + +`centralFinDiffCoeff()` takes no parameters, however it requires the type, order and accuracy as template parameters. + +* `T` - the type you want the coefficients in +* `ord` - the order of the derivative +* `acc` - the accuracy of the approximation + +## Output + +`centralFinDiffCoeff` returns a vector of coefficients. + +## Code +{% highlight c++ %} +template +auto centralFinDiffCoeff() +{ + constexpr int size = 2.0 * std::floor((ord + 1.0) / 2.0) - 1.0 + acc; + constexpr int P = (size - 1.0) / 2.0; + + Matrix mat; + for (auto i = 0; i < size; ++i) + { + for (auto j = 0; j < size ; ++j) + { + mat[i][j] = std::pow(-P+j, i); + } + } + + std::array b; + b.fill(0.0); + b[ord] = fact(ord); + + return mat.solveLinearSystemLU(b); +} +{% endhighlight %} + +## Example +{% highlight c++ %} +int main() +{ + auto coeffs = centralFinDiffCoeff(); + + std::cout << "coefficients of a second order derivative with 4th accuracy\n\n"; + std::cout << coeffs << std::endl; + + return EXIT_SUCCESS; +} +{% endhighlight %} + +## Result +``` +coefficients of a second order derivative with 4th accuracy + +[ -0.0833 1.33 -2.5 1.33 -0.0833 ] + +``` + +**Last Modification date:** 11 January 2018 diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/finiteDiffMethods/manual_init_elliptic_ode.md b/MATH5620_NumericalSolutionsOfDifferentialEquations/finiteDiffMethods/manual_init_elliptic_ode.md new file mode 100644 index 0000000..fbc57b9 --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/finiteDiffMethods/manual_init_elliptic_ode.md @@ -0,0 +1,98 @@ +--- +title: Init Elliptic ODE +math: true +layout: default +--- + +{% include mathjax.html %} + + Table of Contents +# Initialize Elliptic ODE + +**Routine Name:** initEllipticODE + +**Author:** Philip Nelson + +**Language:** C++ + +## Description + +`initEllipticODE` sets up the vectors to solve the elliptic ODE + +## Input + +`initEllipticODE(F f, T a, T b, T ua, T ub)` requires: + +* `F f` - the function \\(f\\) +* `T a` - the left boundary +* `T b` - the right boundary +* `T ua` - the value of \\(u(a)\\) +* `T ub` - the value of \\(u(a)\\) + +## Output + +`initEllipticODE` returns the `a b c` vectors for Thomas Algorithm, which are the diagonals, and the `f`, the function sampled along the mesh. + +## Code +{% highlight c++ %} +template +std::tuple, std::array, std::array, std::array> + initEllipticODE(F f, T a, T b, T ua, T ub) +{ + auto h = (b - a) / N; + auto h2 = h * h; + + auto coeff = centralFinDiffCoeff(); + + std::array av, bv, cv, fv; + + av.fill(coeff[0]); + bv.fill(coeff[1]); + cv.fill(coeff[2]); + + for (auto i = 1u; i < N; ++i) + { + fv[i - 1] = h2 * f(a + i * h); + } + + fv[0] -= ua; + fv[fv.size() - 1] -= ub; + + return {av, bv, cv, fv}; +} +{% endhighlight %} + +## Example +{% highlight c++ %} +int main() +{ + constexpr int N = 10; + auto a = 0.0, b = 1.0, ua = 2.5, ub = 5.0; + auto h = (b - a) / N; + + auto[_a, _b, _c, _f] = initEllipticODE([](double x) { return std::sin(M_PI * x); }, a, b, ua, ub); + + std::cout << "a\n" << _a << std::endl; + std::cout << "b\n" << _b << std::endl; + std::cout << "c\n" << _c << std::endl; + std::cout << "f\n" << _f << std::endl; +} +{% endhighlight %} + +## Result +``` +a +[ 1 1 1 1 1 1 1 1 1 ] + +b +[ -2 -2 -2 -2 -2 -2 -2 -2 -2 ] + +c +[ 1 1 1 1 1 1 1 1 1 ] + +f +[ -2.5 0.00588 0.00809 0.00951 0.01 0.00951 0.00809 0.00588 -5 ] + +``` + +**Last Modification date:** 08 February 2018 diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/finiteDiffMethods/manual_solve_elliptic_ode.md b/MATH5620_NumericalSolutionsOfDifferentialEquations/finiteDiffMethods/manual_solve_elliptic_ode.md new file mode 100644 index 0000000..fec0885 --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/finiteDiffMethods/manual_solve_elliptic_ode.md @@ -0,0 +1,95 @@ +--- +title: Solve Elliptic ODE +math: true +layout: default +--- + +{% include mathjax.html %} + + Table of Contents +# Solve Elliptic ODE + +**Routine Name:** solveEllipticODE + +**Author:** Philip Nelson + +**Language:** C++ + +## Description + +`solveEllipticODE` uses the `initEllipticODE` function and solves the elliptic ODE + +## Input + +`solveEllipticODE(F f, T a, T b, T ua, T ub)` requires: + +* `F f` - the function \\(f\\) +* `T a` - the left boundary +* `T b` - the right boundary +* `T ua` - the value of \\(u(a)\\) +* `T ub` - the value of \\(u(a)\\) + +## Output + +`solveEllipticODE` returns an array of approximations of the elliptic ODE at discrete points from \\(a\\) to \\(b\\) + +## Code +{% highlight c++ %} +template +std::array solveEllipticODE(F f, T a, T b, T ua, T ub) +{ + auto[_a, _b, _c, _f] = initEllipticODE(f, a, b, ua, ub); + + return Matrix::triDiagThomas(_a, _b, _c, _f); +} +{% endhighlight %} + +## Example +{% highlight c++ %} +double exact(double x) +{ + return 2.5 + 2.5 * x - 1 / (M_PI * M_PI) * std::sin(M_PI * x); +} + +int main() +{ + constexpr int N = 10; + auto a = 0.0, b = 1.0, ua = 2.5, ub = 5.0; + auto h = (b - a) / N; + + auto calcVal = + solveEllipticODE([](double x) { return std::sin(M_PI * x); }, a, b, ua, ub); + + + std::array realVal; + for (auto i = 1u; i < N; ++i) + { + realVal[i - 1] = exact(a + i * h); + } + + std::cout << "Calculated Values\n" << calcVal << std::endl; + std::cout << "Real Values\n" << realVal << std::endl; + + std::cout << "P1 norm: " << pNorm(realVal - calcVal, 1) << std::endl; + std::cout << "P2 norm: " << pNorm(realVal - calcVal, 2) << std::endl; + std::cout << "inf norm: " << infNorm(realVal - calcVal) << std::endl; + + return EXIT_SUCCESS; +} +{% endhighlight %} + +## Result +``` + Calculated Values +[ 2.7184 2.94 3.1674 3.4028 3.6478 3.9028 4.1674 4.44 4.7184 ] + + Real Values +[ 2.7187 2.9404 3.168 3.4036 3.6487 3.9036 4.168 4.4404 4.7187 ] + +P1 norm: 0.0052875 +P2 norm: 0.0018726 +inf norm: 0.00083746 + +``` + +**Last Modification date:** 08 February 2018 diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/finiteDiffMethods/manual_solve_elliptic_ode_k.md b/MATH5620_NumericalSolutionsOfDifferentialEquations/finiteDiffMethods/manual_solve_elliptic_ode_k.md new file mode 100644 index 0000000..ecdd341 --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/finiteDiffMethods/manual_solve_elliptic_ode_k.md @@ -0,0 +1,96 @@ +--- +title: Solve Elliptic ODE K +math: true +layout: default +--- + +{% include mathjax.html %} + + Table of Contents +# Solve Elliptic ODE with K + +**Routine Name:** solveEllipticODEwithK + +**Author:** Philip Nelson + +**Language:** C++ + +## Description + +`solveEllipticODEwithK` uses the `initEllipticODE` function and solves the elliptic ODE + +## Input + +`solveEllipticODEwithK(F f, std::array k, T a, T b, T ua, T ub)` requires: + +* `F f` - the function \\(f\\) +* `std::array k` - vector for function \\(k\\) +* `T a` - the left boundary +* `T b` - the right boundary +* `T ua` - the value of \\(u(a)\\) +* `T ub` - the value of \\(u(a)\\) + +## Output + +`solveEllipticODEwithK` returns an array of approximations of the elliptic ODE at discrete points from \\(a\\) to \\(b\\) + +## Code +{% highlight c++ %} +template +std::array solveEllipticODEwithK(F f, T a, T b, T ua, T ub) +{ + auto[_a, _b, _c, _f] = initEllipticODEwithK(f, a, b, ua, ub); + + return Matrix::triDiagThomas(_a, _b, _c, _f); +} +{% endhighlight %} + +## Example +{% highlight c++ %} +double exact(double x) +{ + return 2.5 + 2.5 * x - 1 / (M_PI * M_PI) * std::sin(M_PI * x); +} + +int main() +{ + constexpr int N = 10; + auto a = 0.0, b = 1.0, ua = 2.5, ub = 5.0; + auto h = (b - a) / N; + + auto calcVal = + solveEllipticODEwithK([](double x) { return std::sin(M_PI * x); }, a, b, ua, ub); + + + std::array realVal; + for (auto i = 1u; i < N; ++i) + { + realVal[i - 1] = exact(a + i * h); + } + + std::cout << "Calculated Values\n" << calcVal << std::endl; + std::cout << "Real Values\n" << realVal << std::endl; + + std::cout << "P1 norm: " << pNorm(realVal - calcVal, 1) << std::endl; + std::cout << "P2 norm: " << pNorm(realVal - calcVal, 2) << std::endl; + std::cout << "inf norm: " << infNorm(realVal - calcVal) << std::endl; + + return EXIT_SUCCESS; +} +{% endhighlight %} + +## Result +``` + Calculated Values +[ 2.7184 2.94 3.1674 3.4028 3.6478 3.9028 4.1674 4.44 4.7184 ] + + Real Values +[ 2.7187 2.9404 3.168 3.4036 3.6487 3.9036 4.168 4.4404 4.7187 ] + +P1 norm: 0.0052875 +P2 norm: 0.0018726 +inf norm: 0.00083746 + +``` + +**Last Modification date:** 08 February 2018 diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/gaussSidel/Makefile b/MATH5620_NumericalSolutionsOfDifferentialEquations/gaussSidel/Makefile new file mode 100644 index 0000000..3ceaf8c --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/gaussSidel/Makefile @@ -0,0 +1,13 @@ +OBJS = ../matrix/matrix.hpp ../matrix/matrix_util.hpp ../matrix/vector_util.hpp gaussSidel.hpp main.cpp + +CC = g++ +FLAGS = -std=c++17 +DEBUG_FLAGS = -g3 -Og +RELEASE_FLAGS = -g0 -O3 +WARNINGS = -Wall -Wextra -Werror + +release: $(OBJS) + $(CC) $(FLAGS) $(RELEASE_FLAGS) $(WARNINGS) $(OBJS) -o gaussSidel.out + +debug: $(OBJS) + $(CC) $(FLAGS) $(DEBUG_FLAGS) $(WARNINGS) $(OBJS) -o debug.out diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/gaussSidel/gaussSidel.hpp b/MATH5620_NumericalSolutionsOfDifferentialEquations/gaussSidel/gaussSidel.hpp new file mode 100644 index 0000000..1e7cfd9 --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/gaussSidel/gaussSidel.hpp @@ -0,0 +1,44 @@ +#ifndef GAUSS_SIDEL_HPP +#define GAUSS_SIDEL_HPP + +#include "../machineEpsilon/maceps.hpp" +#include "../matrix/matrix.hpp" +#include "../matrix/vector_util.hpp" +#include + +template +std::array gauss_sidel(Matrix& A, + std::array const& b, + unsigned int const& MAX_ITERATIONS = 1000u) +{ + auto ct = 0u; + std::array x, xn; + x.fill(0); + for (auto k = 0u; k < MAX_ITERATIONS; ++k) + { + ++ct; + xn.fill(0); + for (auto i = 0u; i < N; ++i) + { + auto s1 = 0.0, s2 = 0.0; + for (auto j = 0u; j < i; ++j) + { + s1 += A[i][j] * xn[j]; + } + for (auto j = i + 1; j < N; ++j) + { + s2 += A[i][j] * x[j]; + } + xn[i] = (b[i] - s1 - s2) / A[i][i]; + } + if (allclose(x, xn, maceps().maceps)) + { + std::cout << "Gauss Sidel completed in " << ct << " iterations\n"; + return xn; + } + x = xn; + } + return x; +} + +#endif diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/gaussSidel/main.cpp b/MATH5620_NumericalSolutionsOfDifferentialEquations/gaussSidel/main.cpp new file mode 100644 index 0000000..37f734a --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/gaussSidel/main.cpp @@ -0,0 +1,25 @@ +#include +#include "gaussSidel.hpp" + +int main() +{ + Matrix A ({ + {16, 3}, + {7, -11} + }); + std::array b({{11, 13}}); + std::cout << "A\n" << A << "b\n" << b << '\n'; + std::cout << gauss_sidel(A, b, 100u) << "\n"; + std::cout << A.jacobiIteration(b, 100u) << "\n"; + + Matrix A1 ({ + {10, -1, 2, 0}, + {-1, 11, -1, 3}, + {2, -1, 10, -1}, + {0, 3, -1, 8} + }); + std::array b1({{6, 25, -11, 15}}); + std::cout << "A\n" << A1 << "b\n" << b1 << '\n'; + std::cout << gauss_sidel(A1, b1, 100u) << "\n"; + std::cout << A1.jacobiIteration(b1, 100u) << "\n"; +} diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/gaussSidel/manual_gauss_sidel.md b/MATH5620_NumericalSolutionsOfDifferentialEquations/gaussSidel/manual_gauss_sidel.md new file mode 100644 index 0000000..ea5ade8 --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/gaussSidel/manual_gauss_sidel.md @@ -0,0 +1,138 @@ +--- +title: Gauss Sidel +math: true +layout: default +--- + +{% include mathjax.html %} + + Table of Contents +# Gauss Sidel + +**Routine Name:** gauss_sidel + +**Author:** Philip Nelson + +**Language:** C++ + +## Description + +`gauss_sidel` is an iterative method used to solve a linear system of equations \\(Ax=b\\). It is named after the German mathematicians Carl Friedrich Gauss and Philipp Ludwig von Seidel, and is similar to the Jacobi method. Though it can be applied to any matrix with non-zero elements on the diagonals, convergence is only guaranteed if the matrix is either diagonally dominant, or symmetric and positive definite. [1](https://en.wikipedia.org/wiki/Gauss–Seidel_method) + +## Input + +``` +gauss_sidel(Matrix& A, + std::array const& b, + unsigned int const& MAX_ITERATIONS = 1000u) +``` +requires: + +* `Matrix& A` - an `N`x`N` matrix of type `T` +* `std::array const& b` - the b vector +* `unsigned int const& MAX_ITERATIONS` - the maximum iterations (default 1000) + +## Output + +The solution vector \\(x\\). + +## Code +{% highlight c++ %} +template +std::array gauss_sidel(Matrix& A, + std::array const& b, + unsigned int const& MAX_ITERATIONS = 1000u) +{ + std::array x, xn; + x.fill(0); + for (auto k = 0u; k < MAX_ITERATIONS; ++k) + { + xn.fill(0); + std::cout << x << '\n'; + for (auto i = 0u; i < N; ++i) + { + auto s1 = 0.0, s2 = 0.0; + for (auto j = 0u; j < i; ++j) + { + s1 += A[i][j] * xn[j]; + } + for (auto j = i + 1; j < N; ++j) + { + s2 += A[i][j] * x[j]; + } + xn[i] = (b[i] - s1 - s2) / A[i][i]; + } + if (allclose(x, xn, maceps().maceps)) + { + return xn; + } + x = xn; + } + return x; +} +{% endhighlight %} + +## Example +{% highlight c++ %} +int main() +{ +int main() +{ + Matrix A ({ + {16, 3}, + { 7, -11} + }); + std::array b({ + {11, 13} + }); + std::cout << "A\n" << A << "b\n" << b << '\n'; + std::cout << gauss_sidel(A, b, 100u) << "\n"; + std::cout << A.jacobiIteration(b, 100u) << "\n"; + + Matrix A1({ + {10, -1, 2, 0}, + {-1, 11, -1, 3}, + { 2, -1, 10, -1}, + { 0, 3, -1, 8} + }); + std::array b1({ + {6, 25, -11, 15} + }); + std::cout << "A\n" << A1 << "b\n" << b1 << '\n'; + std::cout << gauss_sidel(A1, b1, 100u) << "\n"; + std::cout << A1.jacobiIteration(b1, 100u) << "\n"; +} +} +{% endhighlight %} + +## Result +It is clear to see that as the matrix size increases, the Gauss Sidel method outperforms Jacobi Iteration. +``` +A +| 16 3 | +| 7 -11 | +b +[ 11 13 ] + +Gauss Sidel completed in 18 iterations +[ 0.812 -0.665 ] + +Jacobi Iteration completed in 35 iterations +[ 0.812 -0.665 ] + +A +| 10 -1 2 0 | +| -1 11 -1 3 | +| 2 -1 10 -1 | +| 0 3 -1 8 | +b +[ 6 25 -11 15 ] + +Gauss Sidel completed in 17 iterations +[ 1 2 -1 1 ] + +Jacobi Iteration completed in 43 iterations +[ 1 2 -1 1 ] +``` + +**Last Modification date:** 20 March 2018 diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/heatEquations/Makefile b/MATH5620_NumericalSolutionsOfDifferentialEquations/heatEquations/Makefile new file mode 100644 index 0000000..c8a0e60 --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/heatEquations/Makefile @@ -0,0 +1,13 @@ +OBJS = heat.hpp main.cpp + +CC = g++ +FLAGS = -std=c++17 +DEBUG_FLAGS = -g3 -Og +RELEASE_FLAGS = -g0 -O3 +WARNINGS = -Wall -Wextra -Werror + +release: $(OBJS) + $(CC) $(FLAGS) $(RELEASE_FLAGS) $(WARNINGS) $(OBJS) -o heat.out + +debug: $(OBJS) + $(CC) $(FLAGS) $(DEBUG_FLAGS) $(WARNINGS) $(OBJS) -o debug.out diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/heatEquations/heat.hpp b/MATH5620_NumericalSolutionsOfDifferentialEquations/heatEquations/heat.hpp new file mode 100644 index 0000000..5761c27 --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/heatEquations/heat.hpp @@ -0,0 +1,54 @@ +#ifndef HEAT_HPP_HPP +#define HEAT_HPP_HPP + +#include +#include + +template +std::vector linspace(T start, T end, unsigned int n) +{ + std::vector x; + auto dx = (end - start) / n; + for (auto i = 0u; i < n; ++i) + { + x.push_back(start + i * dx); + } + + return x; +} + +template +auto heat_explicit_euler(T L, unsigned int Nx, T F, T a, T Tf) +{ + auto x = linspace(0, L, Nx + 1); + auto dx = x[1] - x[0]; + auto dt = F * dx * dx / a; + auto Nt = std::round(Tf / dt); + auto t = linspace(0, Tf, Nt + 1); + std::vector u; + std::vector u_1; + + for (auto i = 0u; i <= Nx; ++i) + { + u.push_back(0); + u_1.push_back(0); + } + + u_1[Nx / 2] = 1000; + + for (auto n = 0u; n < Nt; ++n) + { + for (auto i = 1u; i < Nx; ++i) + { + u[i] = u_1[i] + F * (u_1[i - 1] - 2 * u_1[i] + u_1[i + 1]); + } + + u[0] = 0; + u[Nx] = 0; + std::swap(u, u_1); + } + + return u; +} + +#endif diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/heatEquations/main.cpp b/MATH5620_NumericalSolutionsOfDifferentialEquations/heatEquations/main.cpp new file mode 100644 index 0000000..5407603 --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/heatEquations/main.cpp @@ -0,0 +1,19 @@ +#include "heat.hpp" +#include +#include + +int main() +{ + auto a = 0.7; + auto F = 0.01; + auto Nx = 10u; + auto Tf = 1.0; + auto L = 1.0; + + auto solution = heat_explicit_euler(L, Nx, F, a, Tf); + + for(auto && x:solution) + std::cout << x << '\n'; + + return EXIT_SUCCESS; +} diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/heatEquations/manual_heat_equation.md b/MATH5620_NumericalSolutionsOfDifferentialEquations/heatEquations/manual_heat_equation.md new file mode 100644 index 0000000..3ae707a --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/heatEquations/manual_heat_equation.md @@ -0,0 +1,45 @@ +--- +title: Heat Equation +math: true +layout: default +--- + +{% include mathjax.html %} + + Table of Contents +# Heat Equation + +**Routine Name:** Heat Equation + +**Author:** Philip Nelson + +**Language:** C++ + +## Description + +The heat equation is a linear partial differential equation + +\\[\frac{\partial u}{\partial t} = k \frac{\partial^2 u}{ \partial x^2 }\\] + +that describes the distribution of heat (or variation in temperature) in a given region over time. + +## Input + +## Output + +## Code +{% highlight c++ %} +{% endhighlight %} + +## Example +{% highlight c++ %} +int main() +{ +} +{% endhighlight %} + +## Result +``` +``` + +**Last Modification date:** 5 May 2018 diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/heatEquations/manual_heat_equation_explicit_euler.md b/MATH5620_NumericalSolutionsOfDifferentialEquations/heatEquations/manual_heat_equation_explicit_euler.md new file mode 100644 index 0000000..dda60ac --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/heatEquations/manual_heat_equation_explicit_euler.md @@ -0,0 +1,110 @@ +--- +title: Heat Equations +math: true +layout: default +--- + +{% include mathjax.html %} + + Table of Contents +# Heat Equation + +**Routine Name:** heat_explicit_euler + +**Author:** Philip Nelson + +**Language:** C++ + +## Description + +The heat equation is a linear partial differential equation + +\\[\frac{\partial u}{\partial t} = k \frac{\partial^2 u}{ \partial x^2 }\\] + +that describes the distribution of heat (or variation in temperature) in a given region over time. + +## Input +`heat_explicit_euler(T L, unsigned int Nx, T F, T a, T Tf)` + +* `T L` - Length of the domain ([0,L]) +* `unsigned int Nx` - The total number of mesh points +* `T F` - The dimensionless number \\(\alpha \cdot \frac{dt}{dx^2}\\), which implicitly specifies the time step +* `T a` - Variable coefficient (constant) +* `T Tf` - The stop time for the simulation + +## Output + +A vector containing the state of the system at \\(Tf\\). + +## Code +{% highlight c++ %} +template +auto heat_explicit_euler(T L, unsigned int Nx, T F, T a, T Tf) +{ + auto x = linspace(0, L, Nx + 1); + auto dx = x[1] - x[0]; + auto dt = F * dx * dx / a; + auto Nt = std::round(Tf / dt); + auto t = linspace(0, Tf, Nt + 1); + std::vector u; + std::vector u_1; + + for (auto i = 0u; i <= Nx; ++i) + { + u.push_back(0); + u_1.push_back(0); + } + + u_1[Nx / 2] = 1000; + + for (auto n = 0u; n < Nt; ++n) + { + for (auto i = 1u; i < Nx; ++i) + { + u[i] = u_1[i] + F * (u_1[i - 1] - 2 * u_1[i] + u_1[i + 1]); + } + + u[0] = 0; + u[Nx] = 0; + std::swap(u, u_1); + } + + return u; +} +{% endhighlight %} + +## Example +{% highlight c++ %} +int main() +{ + auto a = 0.7; + auto F = 0.01; + auto Nx = 10u; + auto Tf = 1.0; + auto L = 1.0; + + auto solution = heat_explicit_euler(L, Nx, F, a, Tf); + + for(auto && x:solution) + std::cout << x << '\n'; + + return EXIT_SUCCESS; +} +{% endhighlight %} + +## Result +``` +0 +1.545 +2.93876 +4.04485 +4.75501 +4.99971 +4.75501 +4.04485 +2.93876 +1.545 +0 +``` + +**Last Modification date:** 5 May 2018 diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/heatEquations/manual_heat_equation_implicit_euler.md b/MATH5620_NumericalSolutionsOfDifferentialEquations/heatEquations/manual_heat_equation_implicit_euler.md new file mode 100644 index 0000000..3b61f62 --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/heatEquations/manual_heat_equation_implicit_euler.md @@ -0,0 +1,110 @@ +--- +title: Heat Equations +math: true +layout: default +--- + +{% include mathjax.html %} + + Table of Contents +# Heat Equation Implicit Euler + +**Routine Name:** heat_implicit_euler + +**Author:** Philip Nelson + +**Language:** C++ + +## Description + +The heat equation is a linear partial differential equation + +\\[\frac{\partial u}{\partial t} = k \frac{\partial^2 u}{ \partial x^2 }\\] + +that describes the distribution of heat (or variation in temperature) in a given region over time. + +## Input +`heat_implicit_euler(T L, unsigned int Nx, T F, T a, T Tf)` + +* `T L` - Length of the domain ([0,L]) +* `unsigned int Nx` - The total number of mesh points +* `T F` - The dimensionless number \\(\alpha \cdot \frac{dt}{dx^2}\\), which implicitly specifies the time step +* `T a` - Variable coefficient (constant) +* `T Tf` - The stop time for the simulation + +## Output + +A vector containing the state of the system at \\(Tf\\). + +## Code +{% highlight c++ %} +template +auto heat_implicit_euler(T L, unsigned int Nx, T F, T a, T Tf) +{ + auto x = linspace(0, L, Nx + 1); + auto dx = x[1] - x[0]; + auto dt = F * dx * dx / a; + auto Nt = std::round(Tf / dt); + auto t = linspace(0, Tf, Nt + 1); + std::vector u; + std::vector u_1; + + for (auto i = 0u; i <= Nx; ++i) + { + u.push_back(0); + u_1.push_back(0); + } + + u_1[Nx / 2] = 1000; + + for (auto n = 0u; n < Nt; ++n) + { + for (auto i = 1u; i < Nx; ++i) + { + u[i] = u_1[i] + F * (u_1[i - 1] - 2 * u_1[i] + u_1[i + 1]); + } + + u[0] = 0; + u[Nx] = 0; + std::swap(u, u_1); + } + + return u; +} +{% endhighlight %} + +## Example +{% highlight c++ %} +int main() +{ + auto a = 0.7; + auto F = 0.01; + auto Nx = 10u; + auto Tf = 1.0; + auto L = 1.0; + + auto solution = heat_implicit_euler(L, Nx, F, a, Tf); + + for(auto && x:solution) + std::cout << x << '\n'; + + return EXIT_SUCCESS; +} +{% endhighlight %} + +## Result +``` +0 +1.545 +2.93876 +4.04485 +4.75501 +4.99971 +4.75501 +4.04485 +2.93876 +1.545 +0 +``` + +**Last Modification date:** 5 May 2018 diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/heatEquations/manual_heat_equation_predictor_corrector.md b/MATH5620_NumericalSolutionsOfDifferentialEquations/heatEquations/manual_heat_equation_predictor_corrector.md new file mode 100644 index 0000000..d1e5e33 --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/heatEquations/manual_heat_equation_predictor_corrector.md @@ -0,0 +1,110 @@ +--- +title: Heat Equations +math: true +layout: default +--- + +{% include mathjax.html %} + + Table of Contents +# Heat Equation Predictor Corrector + +**Routine Name:** heat_predictor_corrector + +**Author:** Philip Nelson + +**Language:** C++ + +## Description + +The heat equation is a linear partial differential equation + +\\[\frac{\partial u}{\partial t} = k \frac{\partial^2 u}{ \partial x^2 }\\] + +that describes the distribution of heat (or variation in temperature) in a given region over time. + +## Input +`heat_predictor_corrector(T L, unsigned int Nx, T F, T a, T Tf)` + +* `T L` - Length of the domain ([0,L]) +* `unsigned int Nx` - The total number of mesh points +* `T F` - The dimensionless number \\(\alpha \cdot \frac{dt}{dx^2}\\), which implicitly specifies the time step +* `T a` - Variable coefficient (constant) +* `T Tf` - The stop time for the simulation + +## Output + +A vector containing the state of the system at \\(Tf\\). + +## Code +{% highlight c++ %} +template +auto heat_predictor_corrector(T L, unsigned int Nx, T F, T a, T Tf) +{ + auto x = linspace(0, L, Nx + 1); + auto dx = x[1] - x[0]; + auto dt = F * dx * dx / a; + auto Nt = std::round(Tf / dt); + auto t = linspace(0, Tf, Nt + 1); + std::vector u; + std::vector u_1; + + for (auto i = 0u; i <= Nx; ++i) + { + u.push_back(0); + u_1.push_back(0); + } + + u_1[Nx / 2] = 1000; + + for (auto n = 0u; n < Nt; ++n) + { + for (auto i = 1u; i < Nx; ++i) + { + u[i] = u_1[i] + F * (u_1[i - 1] - 2 * u_1[i] + u_1[i + 1]); + } + + u[0] = 0; + u[Nx] = 0; + std::swap(u, u_1); + } + + return u; +} +{% endhighlight %} + +## Example +{% highlight c++ %} +int main() +{ + auto a = 0.7; + auto F = 0.01; + auto Nx = 10u; + auto Tf = 1.0; + auto L = 1.0; + + auto solution = heat_predictor_corrector(L, Nx, F, a, Tf); + + for(auto && x:solution) + std::cout << x << '\n'; + + return EXIT_SUCCESS; +} +{% endhighlight %} + +## Result +``` +0 +1.545 +2.93876 +4.04485 +4.75501 +4.99971 +4.75501 +4.04485 +2.93876 +1.545 +0 +``` + +**Last Modification date:** 5 May 2018 diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/heatEquations/manual_heat_equation_runge_kutta.md b/MATH5620_NumericalSolutionsOfDifferentialEquations/heatEquations/manual_heat_equation_runge_kutta.md new file mode 100644 index 0000000..1487821 --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/heatEquations/manual_heat_equation_runge_kutta.md @@ -0,0 +1,110 @@ +--- +title: Heat Equation +math: true +layout: default +--- + +{% include mathjax.html %} + + Table of Contents +# Heat Equation + +**Routine Name:** heat_runge_kutta + +**Author:** Philip Nelson + +**Language:** C++ + +## Description + +The heat equation is a linear partial differential equation + +\\[\frac{\partial u}{\partial t} = k \frac{\partial^2 u}{ \partial x^2 }\\] + +that describes the distribution of heat (or variation in temperature) in a given region over time. + +## Input +`heat_runge_kutta(T L, unsigned int Nx, T F, T a, T Tf)` + +* `T L` - Length of the domain ([0,L]) +* `unsigned int Nx` - The total number of mesh points +* `T F` - The dimensionless number \\(\alpha \cdot \frac{dt}{dx^2}\\), which implicitly specifies the time step +* `T a` - Variable coefficient (constant) +* `T Tf` - The stop time for the simulation + +## Output + +A vector containing the state of the system at \\(Tf\\). + +## Code +{% highlight c++ %} +template +auto heat_runge_kutta(T L, unsigned int Nx, T F, T a, T Tf) +{ + auto x = linspace(0, L, Nx + 1); + auto dx = x[1] - x[0]; + auto dt = F * dx * dx / a; + auto Nt = std::round(Tf / dt); + auto t = linspace(0, Tf, Nt + 1); + std::vector u; + std::vector u_1; + + for (auto i = 0u; i <= Nx; ++i) + { + u.push_back(0); + u_1.push_back(0); + } + + u_1[Nx / 2] = 1000; + + for (auto n = 0u; n < Nt; ++n) + { + for (auto i = 1u; i < Nx; ++i) + { + u[i] = u_1[i] + F * (u_1[i - 1] - 2 * u_1[i] + u_1[i + 1]); + } + + u[0] = 0; + u[Nx] = 0; + std::swap(u, u_1); + } + + return u; +} +{% endhighlight %} + +## Example +{% highlight c++ %} +int main() +{ + auto a = 0.7; + auto F = 0.01; + auto Nx = 10u; + auto Tf = 1.0; + auto L = 1.0; + + auto solution = heat_runge_kutta(L, Nx, F, a, Tf); + + for(auto && x:solution) + std::cout << x << '\n'; + + return EXIT_SUCCESS; +} +{% endhighlight %} + +## Result +``` +0 +1.545 +2.93876 +4.04485 +4.75501 +4.99971 +4.75501 +4.04485 +2.93876 +1.545 +0 +``` + +**Last Modification date:** 5 May 2018 diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/hw6_experiments.md b/MATH5620_NumericalSolutionsOfDifferentialEquations/hw6_experiments.md new file mode 100644 index 0000000..d131335 --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/hw6_experiments.md @@ -0,0 +1,99 @@ +--- +title: Homework 6 +math: true +layout: default +--- + +{% include mathjax.html %} + + Table of Contents +# Homework 6 + +**Language:** C++ + +## Description + +Test problems from 7.1, 7.2, and 7.3 of the text. + +{% highlight c++ %} +int main() { + cout << "Problem 1: Explicit Euler" << endl; + cout << Euler(0, 1, 2, .001, .001, 2000, [=](T y, T yk){ return y - yk - h * y; }) << endl; // 7.1 + cout << Euler(0, 1, 2, .001, .001, 2000, [=](T y, T yk){ return y - yk + h * y; }) << endl; // 7.2 + cout << Euler(0, 1, 2, .001, .001, 2000, [=](T y, T yk){ return y - yk + 100 * h * y; }) << endl << endl; // 7.3 + + + cout << "Problem 2: Implicit Euler" << endl; + cout << impEuler(0, 1, 2, .001, .001, 2000, + [=](T y, T yk){ return y-yk - h*y; }, [=](T y, T yk){ return 1 - h; }) << endl; // 7.1 + cout << impEuler(0, 1, 2, .001, .001, 2000, + [=](T y, T yk){ return y- yk + h*y; }, [=](T y, T yk){ return 1 + h; }) << endl; // 7.2 + cout << impEuler(0, 1, 2, .001, .001, 2000, + [=](T y, T yk){ return y- yk + 100*h*y; }, [=](T y, T yk){ return 1 + 100 * h; }) << endl << endl; // 7.3 + + + cout << "Problem 3: Explicit Euler 10^-2" << endl; + cout << Euler(0, 1, 2, .01, .001, 2000, [=](T y, T yk){ return y - yk - h * y; }) << endl; // 7.1 + cout << Euler(0, 1, 2, .01, .001, 2000, [=](T y, T yk){ return y - yk + h * y; }) << endl; // 7.2 + cout << Euler(0, 1, 2, .01, .001, 2000, [=](T y, T yk){ return y - yk + 100 * h * y; }) << endl << endl; // 7.3 + + + cout << "Explicit 10^-1" << endl; + cout << Euler(0, 1, 2, .1, .001, 2000, [=](T y, T yk){ return y - yk - h * y; }) << endl; // 7.1 + cout << Euler(0, 1, 2, .1, .001, 2000, [=](T y, T yk){ return y - yk + h * y; }) << endl; // 7.2 + cout << Euler(0, 1, 2, .1, .001, 2000, [=](T y, T yk){ return y - yk + 100 * h * y; }) << endl << endl; // 7.3 + + + cout << "Problem 4: Adam`s Bashforth - Adam's Moulton" << endl; + cout << AdamBashMoul3(0, 1, 2, .001, [=](T y, T yk){ return y - yk - h * y; }) << endl; // 7.1 + cout << AdamBashMoul3(0, 1, 2, .001, [=](T y, T yk){ return y - yk + h * y; }) << endl; // 7.2 + cout << AdamBashMoul3(0, 1, 2, .001, [=](T y, T yk){ return y - yk + 100 * h * y; }) << endl << endl; // 7.3 + + + cout << "Problem 5: Runge Kutta" << endl; + cout << RungeKutta4(0, 1, 2, .001, [=](T y, T yk){ return y - yk - h * y; }) << endl; // 7.1 + cout << RungeKutta4(0, 1, 2, .001, [=](T y, T yk){ return y - yk + h * y; }) << endl; // 7.2 + cout << RungeKutta4(0, 1, 2, .001, [=](T y, T yk){ return y - yk + 100 * h * y; }) << endl << endl; // 7.3 +} +{% endhighlight %} + +## Results: + +``` +Problem 1: Explicit Euler +-0.415692 +-0.416163 +-1.45252e+76 + + +Problem 2: Implicit Euler +-0.416601 +-0.420292 +-0.833246 + + +Problem 3: Explicit Euler 10^-2 +-0.411589 +-0.416309 +-3.82603e+254 + + +Explicit 10^-1 +-0.369502 +-0.41792 +-6.01633e+41 + + +Problem 4: Adam`s Bashforth - Adam's Moulton +-0.416147 +-0.416147 +7.69951e+283 + + +Problem 5: Runge Kutta +-0.416146 +-0.416147 +-0.416147 +``` + +**Last Modification date:** 3 April 2018 diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/implicitEuler/Makefile b/MATH5620_NumericalSolutionsOfDifferentialEquations/implicitEuler/Makefile new file mode 100644 index 0000000..c9e7068 --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/implicitEuler/Makefile @@ -0,0 +1,13 @@ +OBJS = ../matrix/matrix.hpp ../matrix/matrix_util.hpp ../matrix/vector_util.hpp main.cpp + +CC = g++ +FLAGS = -std=c++17 +DEBUG_FLAGS = -g3 -Og +RELEASE_FLAGS = -g0 -O3 +WARNINGS = -Wall -Wextra -Werror + +release: $(OBJS) + $(CC) $(FLAGS) $(RELEASE_FLAGS) $(WARNINGS) $(OBJS) -o implicitEuler.out + +debug: $(OBJS) + $(CC) $(FLAGS) $(DEBUG_FLAGS) $(WARNINGS) $(OBJS) -o debug.out diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/implicitEuler/implicitEuler.hpp b/MATH5620_NumericalSolutionsOfDifferentialEquations/implicitEuler/implicitEuler.hpp new file mode 100644 index 0000000..b2a547d --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/implicitEuler/implicitEuler.hpp @@ -0,0 +1,25 @@ +#ifndef IMPLICIT_EULER_HPP +#define IMPLICIT_EULER_HPP + +#include "../machineEpsilon/maceps.hpp" +#include "../newtonsMethod/newtonsMethod.hpp" +#include +#include + +template +implicit_euler (F f,T df,T x0,T t0,T tf,const unsigned int MAX_ITERATIONS) +{ + auto h=(tf-t0)/n; + auto t=linspace(t0,tf,n+1); + auto y=zeros(n+1,length(y0)); + auto y(1,:)=y0; + for (auto i=0u; i < MAX_ITERATIONS; ++i) + { + x0=y(i,:)’; + x1=x0-inv(eye(length(y0))-h*feval(df,t(i),x0))*(x0-h*feval(f,t(i),x0)’-y(i,:)’); + x1 = newtons_method(f, dt, x0) + y(i+1,:)=x1’; + } +} + +#endif diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/implicitEuler/main.cpp b/MATH5620_NumericalSolutionsOfDifferentialEquations/implicitEuler/main.cpp new file mode 100644 index 0000000..e97758e --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/implicitEuler/main.cpp @@ -0,0 +1,6 @@ +#include "explicitEuler.hpp" +#include + +int main() { + std::cout << explicit_euler(0.0, -1.0, 0.4, 0.1, [](double a, double b){return a*a+2*b;}) << std::endl; +} diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/implicitEuler/manual_implicit_euler.md b/MATH5620_NumericalSolutionsOfDifferentialEquations/implicitEuler/manual_implicit_euler.md new file mode 100644 index 0000000..ff221c3 --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/implicitEuler/manual_implicit_euler.md @@ -0,0 +1,107 @@ +--- +title: Implicit Euler +math: true +layout: default +--- + +{% include mathjax.html %} + + Table of Contents +# Implicit Euler Method + +**Routine Name:** `implicit_euler` + +**Author:** Philip Nelson + +**Language:** C++ + +## Description + +`implicit_euler`, also known as the backward Euler method, is one of the most basic numerical methods for the solution of ordinary differential equations. It is similar to the (standard) Euler method, but differs in that it is an implicit method. The backward Euler method has order one in time. [1](https://en.wikipedia.org/wiki/Backward_Euler_method) + +## Code +{% highlight c++ %} +template +implicit_euler (F f,T df,T x0,T t0,T tf,const unsigned int MAX_ITERATIONS) +{ + auto h=(tf-t0)/n; + auto t=linspace(t0,tf,n+1); + auto y=zeros(n+1,length(y0)); + auto y(1,:)=y0; + for (auto i=0u; i < MAX_ITERATIONS; ++i) + { + x0=y(i,:); + x1=x0-inv(eye(length(y0))-h*feval(df,t(i),x0))*(x0-h*feval(f,t(i),x0)’-y(i,:)’); + x1 = newtons_method(f, dt, x0) + y(i+1,:)=x1’; + } +} +{% endhighlight %} + +## Example +{% highlight c++ %} +int main() +{ + std::cout << + implicit_euler(0.0, -1.0, 0.4, 0.1, [](double a, double b){return a*a+2*b;}) + << std::endl; +} +{% endhighlight %} + +## Result +``` +----- Lambda Differential Equation ----- + +lambda = 1 +exact 0 -> 10 +approx 0 -> 10 +exact 0.2 -> 12.214 +approx 0.2 -> 12.214 +exact 0.4 -> 14.9182 +approx 0.4 -> 14.9183 +exact 0.6 -> 18.2212 +approx 0.6 -> 18.2212 + +lambda = -1 +exact 0 -> 10 +approx 0 -> 10 +exact 0.2 -> 8.18731 +approx 0.2 -> 8.18732 +exact 0.4 -> 6.7032 +approx 0.4 -> 6.70321 +exact 0.6 -> 5.48812 +approx 0.6 -> 5.48813 + +lambda = 100 +exact 0 -> 10 +approx 0 -> 10 +exact 0.2 -> 4.85165e+09 +approx 0.2 -> 4.90044e+09 +exact 0.4 -> 2.35385e+18 +approx 0.4 -> inf +exact 0.6 -> 1.14201e+27 + +----- Logistic Differential Equaiton ----- + +p0 = 25 +exact 0 -> 25 +approx 0 -> 25 +exact 0.2 -> 25.4922 +approx 0.2 -> 25.4922 +exact 0.4 -> 25.9937 +approx 0.4 -> 25.9937 +exact 0.6 -> 26.5049 +approx 0.6 -> 26.5049 + +p0 = 40000 +exact 0 -> 40000 +approx 0 -> 40000 +exact 0.2 -> 22570.2 +approx 0.2 -> 22570.4 +exact 0.4 -> 15815.2 +approx 0.4 -> 15815.4 +exact 0.6 -> 12228 +approx 0.6 -> 12228.2 +``` + +**Last Modification date:** 3 April 2018 diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/jacobiIteration/Makefile b/MATH5620_NumericalSolutionsOfDifferentialEquations/jacobiIteration/Makefile new file mode 100644 index 0000000..92d4156 --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/jacobiIteration/Makefile @@ -0,0 +1,12 @@ +OBJS = main.cpp + +CC = g++ +FLAGS = -std=c++17 +DEBUG_FLAGS = -Og -g3 -Wall -Wextra +RELEASE_FLAGS = -O3 + +release: $(OBJS) + $(CC) $(RELEASE_FLAGS) $(FLAGS) $(OBJS) -o jacobiIteration.out + +debug: $(OBJS) + $(CC) $(DEBUG_FLAGS) $(FLAGS) $(OBJS) -o debug.out diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/jacobiIteration/jacobiIteration.hpp b/MATH5620_NumericalSolutionsOfDifferentialEquations/jacobiIteration/jacobiIteration.hpp new file mode 100644 index 0000000..cc1c255 --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/jacobiIteration/jacobiIteration.hpp @@ -0,0 +1,47 @@ +#ifndef JACOBI_ITERATION_HPP +#define JACOBI_ITERATION_HPP + +#include "../matrix/matrix.hpp" +#include "../matrix/vector_util.hpp" +#include + +template +std::array jacobiIteration(Matrix A, + std::array const& b, + unsigned int const& MAX_ITERATIONS = 1000u) +{ + auto ct = 0u; + std::array zeros; + zeros.fill(0); + + std::array x = zeros; + + for (auto n = 0u; n < MAX_ITERATIONS; ++n) + { + ++ct; + auto x_n = zeros; + + for (auto i = 0u; i < N; ++i) + { + T sum = 0; + for (auto j = 0u; j < N; ++j) + { + if (j == i) continue; + sum += A[i][j] * x[j]; + } + x_n[i] = (b[i] - sum) / A[i][i]; + } + + if (allclose(x, x_n, maceps().maceps)) + { + std::cout << "Jacobi Iteration completed in " << ct << " iterations\n"; + return x_n; + } + + x = x_n; + } + + return x; +} + +#endif diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/jacobiIteration/main.cpp b/MATH5620_NumericalSolutionsOfDifferentialEquations/jacobiIteration/main.cpp new file mode 100644 index 0000000..962ae6e --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/jacobiIteration/main.cpp @@ -0,0 +1,21 @@ +#include "../matrix/matrix.hpp" +#include "jacobiIteration.hpp" +#include + +int main() +{ + Matrix A({ + {-2, 1, 0, 0}, + { 1, -2, 1, 0}, + { 0, 1, -2, 1}, + { 0, 0, 1, -2} + }); + std::array x = { + {4, 7, 2, 5} + }; + auto b = A * x; + std::cout << " A\n" << A << std::endl; + std::cout << " b\n" << b << std::endl; + std::cout << " x\n" << x << std::endl; + std::cout << jacobiIteration(A, b, 1000u) << std::endl; +} diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/jacobiIteration/manual_jacobi_iteration.md b/MATH5620_NumericalSolutionsOfDifferentialEquations/jacobiIteration/manual_jacobi_iteration.md new file mode 100644 index 0000000..b64ec4c --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/jacobiIteration/manual_jacobi_iteration.md @@ -0,0 +1,121 @@ +--- +title: Jacobi Iteration +math: true +layout: default +--- + +{% include mathjax.html %} + + Table of Contents +# Jacobi Iteration + +**Routine Name:** jacobiIteration + +**Author:** Philip Nelson + +**Language:** C++ + +## Description + +`jacobiIteration` solves a linear system of equations by Jacobi Iteration. + +## Input + +``` +std::array jacobiIteration(Matrix A, + std::array const& b, + unsigned int const& MAX = 1000u) +``` +requires: + +* `Matrix& A` - an `N`x`N` matrix of type `T` +* `std::array const& b` - the b vector +* `unsigned int const& MAX_ITERATIONS` - the maximum iterations (default 1000) + +## Output + +`jacobiIteration` returns a `std::array` with the solution vector `x` + +## Code +{% highlight c++ %} +template +std::array jacobiIteration(Matrix A, + std::array const& b, + unsigned int const& MAX = 1000u) +{ + auto ct = 0u; + std::array zeros; + zeros.fill(0); + + std::array x = zeros; + + for (auto n = 0u; n < MAX; ++n) + { + ++ct; + auto x_n = zeros; + + for (auto i = 0u; i < N; ++i) + { + T sum = 0; + for (auto j = 0u; j < N; ++j) + { + if (j == i) continue; + sum += A[i][j] * x[j]; + } + x_n[i] = (b[i] - sum) / A[i][i]; + } + + if (allclose(x, x_n, maceps().maceps)) + { + std::cout << "Jacobi Iteration completed in " << ct << " iterations\n"; + return x_n; + } + + x = x_n; + } + + return x; +} +{% endhighlight %} + +## Example +{% highlight c++ %} +int main() +{ + Matrix A({ + {-2, 1, 0, 0}, + { 1, -2, 1, 0}, + { 0, 1, -2, 1}, + { 0, 0, 1, -2} + }); + std::array x = { + {4, 7, 2, 5} + }; + + auto b = A * x; + std::cout << " A\n" << A << std::endl; + std::cout << " b\n" << b << std::endl; + std::cout << " x\n" << x << std::endl; + std::cout << jacobiIteration(A, b, 1000u) << std::endl; +} +{% endhighlight %} + +## Result +``` + A +| -2 1 0 0 | +| 1 -2 1 0 | +| 0 1 -2 1 | +| 0 0 1 -2 | + + b +[ -1 -8 8 -8 ] + + x +[ 4 7 2 5 ] + +Jacobi Iteration completed in 170 iterations +[ 4 7 2 5 ] +``` + +**Last Modification date:** 07 February 2018 diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/laxWendroff/Makefile b/MATH5620_NumericalSolutionsOfDifferentialEquations/laxWendroff/Makefile new file mode 100644 index 0000000..c5d3f59 --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/laxWendroff/Makefile @@ -0,0 +1,13 @@ +OBJS = ../matrix/matrix.hpp ../matrix/matrix_util.hpp ../matrix/vector_util.hpp main.cpp + +CC = g++ +FLAGS = -std=c++17 +DEBUG_FLAGS = -g3 -Og +RELEASE_FLAGS = -g0 -O3 +WARNINGS = -Wall -Wextra -Werror + +release: $(OBJS) + $(CC) $(FLAGS) $(RELEASE_FLAGS) $(WARNINGS) $(OBJS) -o laxWendroff.out + +debug: $(OBJS) + $(CC) $(FLAGS) $(DEBUG_FLAGS) $(WARNINGS) $(OBJS) -o debug.out diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/laxWendroff/laxWendroff.hpp b/MATH5620_NumericalSolutionsOfDifferentialEquations/laxWendroff/laxWendroff.hpp new file mode 100644 index 0000000..96c5733 --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/laxWendroff/laxWendroff.hpp @@ -0,0 +1,42 @@ +#ifndef LAX_WENDORF_HPP +#define LAX_WENDORF_HPP + +#include "../matrix/matrix.hpp" +#include "../matrix/matrix_util.hpp" +#include +#include + +template +std::vector> lax_wendorf(const T xDomain[], + const T tDomain[], + const T dx, + const T dt, + F eta, + const T c) +{ + std::array U; + for (auto i = 0u; i < S - 2; ++i) + { + auto A = xDomain[0]; + auto x = A + (i + 1) * dx; + U[i] = eta(x); + } + + auto a = c * dx / dt / dt; + auto b = c * c * dx * dx / dt / dt / 2; + std::vector coeffs = {a + b, 1 - 2 * b, -a + b}; + + auto A = makeNDiag(coeffs, -1u); + + std::vector> solution = {U}; + + for (auto t = tDomain[0]; t < tDomain[1]; t += dt) + { + U = A * U; + solution.push_back(U); + } + + return solution; +} + +#endif diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/laxWendroff/main.cpp b/MATH5620_NumericalSolutionsOfDifferentialEquations/laxWendroff/main.cpp new file mode 100644 index 0000000..fa471cd --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/laxWendroff/main.cpp @@ -0,0 +1,32 @@ +#include "laxWendroff.hpp" +#include +#include + +int main() +{ + constexpr double xDomain[] = {0.0, 1.0}; + constexpr double tDomain[] = {0.0, 1.0}; + // constexpr auto dt = 1.0e-3; + // constexpr auto dx = 1.0e-3; + constexpr auto dt = .1; + constexpr auto dx = .1; + constexpr std::size_t S = ((xDomain[1] - xDomain[0]) / dx); + auto c = 0.7; + auto eta = [](const double& x) -> double { + return (x >= 0.3 && x <= 0.6) ? 100 : 0; + }; + + auto solution = lax_wendorf(xDomain, tDomain, dx, dt, eta, c); + + for (auto i = 0u; i < solution.size(); ++i) + { + for (auto j = 0u; j < solution[i].size(); ++j) + { + std::cout << std::setprecision(3) << std::setw(10) << solution[j][i] + << " "; + } + std::cout << '\n'; + } + + return EXIT_SUCCESS; +} diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/laxWendroff/manual_lax_wendroff.md b/MATH5620_NumericalSolutionsOfDifferentialEquations/laxWendroff/manual_lax_wendroff.md new file mode 100644 index 0000000..8d451e4 --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/laxWendroff/manual_lax_wendroff.md @@ -0,0 +1,129 @@ +--- +title: Lax Wendroff +math: true +layout: default +--- + +{% include mathjax.html %} + + Table of Contents +# Lax Wendroff + +**Routine Name:** `lax_wendorf` + +**Author:** Philip Nelson + +**Language:** C++ + +## Description + +`lax_wendorf` is a numerical method for the solution of hyperbolic partial differential equations, based on finite differences. It is second-order accurate in both space and time. This method is an example of explicit time integration where the function that defines governing equation is evaluated at the current time.[1](https://en.wikipedia.org/wiki/Lax–Wendroff_method) + +## Input + +{% highlight c++ %} +lax_wendorf(const T xDomain[], + const T tDomain[], + const T dx, + const T dt, + F eta, + const T c + ) +{% endhighlight %} + +* `T xDomain[]` - two member array with the spacial bounds +* `T tDomain[]` - two member array with the temporal bounds +* `T dx` - the spacial step +* `T dt` - the temporal step +* `F eta` - the function eta +* `T c` - the constant + +## Output + +A `std::vector>` with the solution over time. + +## Code +{% highlight c++ %} +template +std::vector> lax_wendorf(const T xDomain[], + const T tDomain[], + const T dx, + const T dt, + F eta, + const T c) +{ + std::array U; + for (auto i = 0u; i < S - 2; ++i) + { + auto A = xDomain[0]; + auto x = A + (i + 1) * dx; + U[i] = eta(x); + } + + auto a = c * dx / dt / dt; + auto b = c * c * dx * dx / dt / dt / 2; + std::vector coeffs = {a + b, 1 - 2 * b, -a + b}; + + auto A = makeNDiag(coeffs, -1u); + + std::vector> solution = {U}; + + for (auto t = tDomain[0]; t < tDomain[1]; t += dt) + { + U = A * U; + solution.push_back(U); + } + + return solution; +} +{% endhighlight %} + +## Example +{% highlight c++ %} +int main() +{ + constexpr double xDomain[] = {0.0, 1.0}; + constexpr double tDomain[] = {0.0, 1.0}; + // constexpr auto dt = 1.0e-3; + // constexpr auto dx = 1.0e-3; + constexpr auto dt = .1; + constexpr auto dx = .1; + constexpr std::size_t S = ((xDomain[1] - xDomain[0]) / dx); + auto c = 0.7; + auto eta = [](const double& x) -> double { + return (x >= 0.3 && x <= 0.6) ? 100 : 0; + }; + + auto solution = lax_wendorf(xDomain, tDomain, dx, dt, eta, c); + + for (auto i = 0u; i < solution.size(); ++i) + { + for (auto j = 0u; j < solution[i].size(); ++j) + { + std::cout << std::setprecision(3) << std::setw(10) << solution[j][i] + << " "; + } + std::cout << '\n'; + } + + return EXIT_SUCCESS; +} +{% endhighlight %} + +## Result +``` + 0 0 4.56e+03 -2.38e+04 -5.17e+05 4.78e+06 6.4e+07 -8.34e+08 + 0 -675 3.87e+03 7.48e+04 -7.47e+05 -9.11e+06 1.28e+08 1.24e+09 + 100 -624 -5.89e+03 9.07e+04 7.38e+05 -1.45e+07 -1.06e+08 2.4e+09 + 100 100 -9.71e+03 -2.21e+04 1.41e+06 4.81e+06 -2.25e+08 -9.4e+08 + 100 775 -3.77e+03 -1.13e+05 1.86e+05 1.81e+07 8.56e+06 -3.04e+09 + 0 724 5.99e+03 -5.97e+04 -1.16e+06 5.26e+06 2.09e+08 -2.62e+08 + 0 0 5.25e+03 4.61e+04 -6.66e+05 -1.11e+07 6.38e+07 2.11e+09 + 0 0 0 3.8e+04 3.53e+05 -4.65e+06 -8.29e+07 4.2e+08 + 0 4.56e+03 -2.38e+04 -5.17e+05 4.78e+06 6.4e+07 -8.34e+08 -8.83e+09 + -675 3.87e+03 7.48e+04 -7.47e+05 -9.11e+06 1.28e+08 1.24e+09 -2.16e+10 + -624 -5.89e+03 9.07e+04 7.38e+05 -1.45e+07 -1.06e+08 2.4e+09 1.66e+10 + 100 -9.71e+03 -2.21e+04 1.41e+06 4.81e+06 -2.25e+08 -9.4e+08 3.74e+10 +``` + +**Last Modification date:** 5 May 2018 diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/logistic/manual.md b/MATH5620_NumericalSolutionsOfDifferentialEquations/logistic/manual.md index eb93979..468b075 100644 --- a/MATH5620_NumericalSolutionsOfDifferentialEquations/logistic/manual.md +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/logistic/manual.md @@ -6,8 +6,8 @@ layout: default {% include mathjax.html %} - Table of Contents -# Logistic Differential Equaiton + Table of Contents +# Logistic Differential Equation **Routine Name:** logistic diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/logistic2/Makefile b/MATH5620_NumericalSolutionsOfDifferentialEquations/logistic2/Makefile new file mode 100644 index 0000000..5504db8 --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/logistic2/Makefile @@ -0,0 +1,12 @@ +OBJS = main.cpp logisticSolver.hpp + +CC = g++ +FLAGS = -std=c++17 +DEBUG_FLAGS = -Og -g3 -Wall -Wextra -Werror +RELEASE_FLAGS = -O3 + +release: $(OBJS) + $(CC) $(RELEASE_FLAGS) $(FLAGS) $(OBJS) -o logistic.out + +debug: $(OBJS) + $(CC) $(DEBUG_FLAGS) $(FLAGS) $(OBJS) -o debug.out diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/logistic2/logisticSolver.hpp b/MATH5620_NumericalSolutionsOfDifferentialEquations/logistic2/logisticSolver.hpp new file mode 100644 index 0000000..ad8b9df --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/logistic2/logisticSolver.hpp @@ -0,0 +1,14 @@ +#ifndef LOGISTIC_SOLVER_HPP +#define LOGISTIC_SOLVER_HPP + +#include + +template +auto logisticSolver(T const& a, T const& b, T const& p0) +{ + return [=](T const& t) { + return (a / (((a - p0 * b) / p0) * std::exp(-a * t) + b)); + }; +} + +#endif diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/logistic2/main.cpp b/MATH5620_NumericalSolutionsOfDifferentialEquations/logistic2/main.cpp new file mode 100644 index 0000000..d4ec595 --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/logistic2/main.cpp @@ -0,0 +1,21 @@ +#include "logisticSolver.hpp" +#include + +int main() +{ + double a = 2.5; + double b = 1.7; + double p0 = 3.2; + + std::cout << "alpha:\t" << a << "\nbeta:\t" << b << "\nP0:\t" << p0 << "\n\n"; + + auto solveLog = logisticSolver(a, b, p0); + + // Call it for some basic values + for (int t = -10; t < 0; ++t) { + std::cout << "t = " << t << " -> " << solveLog(t) + << "\tt = " << t+10 << " -> " << solveLog(t+10) << '\n'; + } + + return EXIT_SUCCESS; +} diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/logistic2/manual.md b/MATH5620_NumericalSolutionsOfDifferentialEquations/logistic2/manual.md new file mode 100644 index 0000000..854c2cf --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/logistic2/manual.md @@ -0,0 +1,90 @@ +--- +title: LogisticSolver +math: true +layout: default +--- + +{% include mathjax.html %} + + Table of Contents +# Logistic Differential Equation + +**Routine Name:** logisticSolver + +**Author:** Philip Nelson + +**Language:** C++ + +## Description + +`logisticSolver` generates a function that takes time and returns the solution to the exact solution to the logistic differential equation: + +\\[ \frac{dP}{dt} = \alpha P + \beta P^2 \\] + +The exact solution has the form: + +\\[\frac{\alpha}{(\frac{\alpha}{P_0}-\Beta)e^{-\alpha t}+\Beta}\\] + +## Input + +`logisticSolver(T const& a, T const& b, T const& p0)` requires: +* `a` - \\(\alpha\\) +* `b` - \\(\beta\\) +* `p0` - initial condition \\(P(0)\\) + +## Output + +Logistic returns a function that takes a \\(t\\) and evaluates the logistic equation at that time. + +## Code +{% highlight c++ %} +template +auto logisticSolver(T const& a, T const& b, T const& p0) +{ + return [=](T const& t) { + return (a / (((a - p0 * b) / p0) * std::exp(-a * t) + b)); + }; +} +{% endhighlight %} + +## Example +{% highlight c++ %} +int main() +{ + double a = 2.5; + double b = 1.7; + double p0 = 3.2; + + std::cout << "alpha:\t" << a << "\nbeta:\t" << b << "\nP0:\t" << p0 << "\n\n"; + + auto solveLog = logisticSolver(a, b, p0); + + // Call it for some basic values + for (int t = -10; t < 0; ++t) { + std::cout << "t = " << t << " -> " << solveLog(t) + << "\tt = " << t+10 << " -> " << solveLog(t+10) << '\n'; + } + + return EXIT_SUCCESS; +} +{% endhighlight %} + +## Result +``` +alpha: 2.5 +beta: 1.7 +P0: 3.2 + +t = -10 -> -3.77903e-11 t = 0 -> 3.2 +t = -9 -> -4.6038e-10 t = 1 -> 1.53886 +t = -8 -> -5.60858e-09 t = 2 -> 1.47596 +t = -7 -> -6.83265e-08 t = 3 -> 1.47103 +t = -6 -> -8.32388e-07 t = 4 -> 1.47062 +t = -5 -> -1.01406e-05 t = 5 -> 1.47059 +t = -4 -> -0.000123548 t = 6 -> 1.47059 +t = -3 -> -0.00150653 t = 7 -> 1.47059 +t = -2 -> -0.018566 t = 8 -> 1.47059 +t = -1 -> -0.263361 t = 9 -> 1.47059 +``` + +**Last Modification date:** 3 April 2018 diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/machineEpsilon/maceps.hpp b/MATH5620_NumericalSolutionsOfDifferentialEquations/machineEpsilon/maceps.hpp index db6cee0..f9dd9d0 100644 --- a/MATH5620_NumericalSolutionsOfDifferentialEquations/machineEpsilon/maceps.hpp +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/machineEpsilon/maceps.hpp @@ -1,3 +1,6 @@ +#ifndef MACEPS_HPP +#define MACEPS_HPP + #include template struct eps @@ -22,3 +25,5 @@ eps maceps() return eps(prec, e); } + +#endif diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/machineEpsilon/manual.md b/MATH5620_NumericalSolutionsOfDifferentialEquations/machineEpsilon/manual.md index 41b10b8..c5608da 100644 --- a/MATH5620_NumericalSolutionsOfDifferentialEquations/machineEpsilon/manual.md +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/machineEpsilon/manual.md @@ -2,7 +2,7 @@ title: Machine Epsilon layout: default --- - Table of Contents + Table of Contents # Machine Epsilon **Routine Name:** maceps diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/manualTemplate.md b/MATH5620_NumericalSolutionsOfDifferentialEquations/manualTemplate.md deleted file mode 100644 index ac2a4d2..0000000 --- a/MATH5620_NumericalSolutionsOfDifferentialEquations/manualTemplate.md +++ /dev/null @@ -1,68 +0,0 @@ - Table of Contents -# Machine Epsilon - -**Routine Name:** maceps - -**Author:** Philip Nelson - -**Language:** C++ - -## Description - -maceps returns the machine epsilon and precision of any primitive type. A make file is included with a driver program. - -``` -$ make -$ ./maceps.out -``` - -This will compile and run the driver program. - -## Input - -`maceps( )` requires a template argument _T_ with the type of machine epsilon you want _( float, double, long double, etc... )_. Otherwise, `maceps` takes no input. - -## Output - -maceps returns an `eps` struct with members `int prec` which holds the precision and `T maceps` which holds the machine epsilon for the specified type. - -## Code -{% highlight c++ %} -template -eps maceps() -{ - T e = 1; - T one = 1; - T half = 0.5; - int prec = 1; - while (one + e * half > one) - { - e *= half; - ++prec; - } - - return eps(prec, e); -} -{% endhighlight %} - -## Example -{% highlight c++ %} -int main() -{ - auto doubleeps = maceps(); - std::cout << "double\n"; - std::cout << "precision:\t" << doubleeps.prec << std::endl; - std::cout << "maceps:\t\t" << doubleeps.maceps << std::endl; - std::cout << "std::numeric:\t" << std::numeric_limits::epsilon() << std::endl << std::endl; -} -{% endhighlight %} - -## Result -``` -double -precision: 53 -maceps: 2.22045e-16 -std::numeric: 2.22045e-16 -``` - -**Last Modification date:** 11 January 2018 diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/Makefile b/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/Makefile index 3b149d1..3c06ac7 100644 --- a/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/Makefile +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/Makefile @@ -2,7 +2,7 @@ OBJS = main.cpp matrix.hpp CC = g++ FLAGS = -std=c++17 -DEBUG_FLAGS = -Og -g3 -Wall -Wextra -Werror +DEBUG_FLAGS = -Og -g3 -Wall -Wextra RELEASE_FLAGS = -O3 release: $(OBJS) diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/example_inverse_power_iteration_elliptic_ode.md b/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/example_inverse_power_iteration_elliptic_ode.md new file mode 100644 index 0000000..c298058 --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/example_inverse_power_iteration_elliptic_ode.md @@ -0,0 +1,76 @@ +--- +title: Inverse Power Iteration Example on the Elliptic ODE +math: true +layout: default +--- + +{% include mathjax.html %} + + Table of Contents +# Example of Inverse Power Iteration + +**Author:** Philip Nelson + +**Language:** C++ + +## Description + +This demonstrates the condition number method on the tri-diagonal matrix used to solve the Elliptic ODE. + +## Output + +The condition number of the `N`x`N` matrix with increasing `N`. Observe that as the size of the matrix increases, the condition number increases unbounded + +## Example +{% highlight c++ %} +int main() +{ + std::cout << "A\n" << SecOrdFinDifMethEllipticMat() << std::endl << std::endl; + std::cout << "3x3" << std::endl; + std::cout << conditionNumber(SecOrdFinDifMethEllipticMat()) << std::endl << std::endl; + std::cout << "5x5" << std::endl; + std::cout << conditionNumber(SecOrdFinDifMethEllipticMat()) << std::endl << std::endl; + std::cout << "10x10" << std::endl; + std::cout << conditionNumber(SecOrdFinDifMethEllipticMat()) << std::endl << std::endl; + std::cout << "15x15" << std::endl; + std::cout << conditionNumber(SecOrdFinDifMethEllipticMat()) << std::endl << std::endl; + std::cout << "25x25" << std::endl; + std::cout << conditionNumber(SecOrdFinDifMethEllipticMat()) << std::endl << std::endl; + std::cout << "75x75" << std::endl; + std::cout << conditionNumber(SecOrdFinDifMethEllipticMat()) << std::endl << std::endl; + std::cout << "100x100" << std::endl; + std::cout << conditionNumber(SecOrdFinDifMethEllipticMat()) << std::endl << std::endl; +} +{% endhighlight %} + +## Result +``` +A +| -2 1 0 | +| 1 -2 1 | +| 0 1 -2 | + +3x3 +-5.82843 + +5x5 +-13.9282 + +10x10 +-48.3742 + +15x15 +-103.087 + +25x25 +-273.306 + +75x75 +-2340.09 + +100x100 +-4133.2 + +``` + +**Last Modification date:** 27 February 2018 diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/example_power_iteration_elliptic_ode.md b/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/example_power_iteration_elliptic_ode.md new file mode 100644 index 0000000..bb07c15 --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/example_power_iteration_elliptic_ode.md @@ -0,0 +1,60 @@ +--- +title: Power Iteration Example on the Elliptic ODE +math: true +layout: default +--- + +{% include mathjax.html %} + + Table of Contents +# Example of Power Iteration + +**Author:** Philip Nelson + +**Language:** C++ + +## Description + +This demonstrates the power iteration method on the tri-diagonal matrix used to solve the Elliptic ODE. + +## Output + +The largest Eigenvalue of the `N`x`N` matrix with increasing `N`. Observe that as the size of the matrix increases, the largest Eigenvalue asymptotically approaches the value 4. + +## Example +{% highlight c++ %} +int main() +{ + std::cout << "A\n" << SecOrdFinDifMethEllipticMat() << std::endl << std::endl; + std::cout << "5x5" << std::endl; + std::cout << powerIteration(SecOrdFinDifMethEllipticMat(), 1000u) << std::endl << std::endl; + std::cout << "10x10" << std::endl; + std::cout << powerIteration(SecOrdFinDifMethEllipticMat(), 1000u) << std::endl << std::endl; + std::cout << "100x100" << std::endl; + std::cout << powerIteration(SecOrdFinDifMethEllipticMat(), 1000u) << std::endl << std::endl; + std::cout << "1000x10000" << std::endl; + std::cout << powerIteration(SecOrdFinDifMethEllipticMat(), 1000u) << std::endl; +} +{% endhighlight %} + +## Result +``` +A +| -2 1 0 | +| 1 -2 1 | +| 0 1 -2 | + +5x5 +3.73205 + +10x10 +3.91899 + +100x100 +3.99852 + +1000x10000 +3.99893 +``` + +**Last Modification date:** 27 February 2018 diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/main.cpp b/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/main.cpp index 225ac37..5b11fe2 100644 --- a/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/main.cpp +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/main.cpp @@ -1,6 +1,10 @@ +#include "../finiteDiffMethods/finDiffCoeff.hpp" +#include "../jacobiIteration/jacobiIteration.hpp" #include "matrix.hpp" #include "matrix_util.hpp" #include "termColors.hpp" +#include "vector_util.hpp" +#include #include #include @@ -14,7 +18,7 @@ void test(T a, R r, std::string name) std::cout << RED << "[ FAIL] " << RESET << std::endl; } -int main() +void runTests() { Matrix a({{4, 0}, {1, -9}}); Matrix _a({{8, 0}, {2, -18}}); @@ -28,14 +32,22 @@ int main() Matrix _e_neg({{-4, 0}, {-1, 9}}); Matrix _id({{1, 0, 0}, {0, 1, 0}, {0, 0, 1}}); Matrix f({{6, 1, 1}, {4, -2, 5}, {2, 8, 7}}); - Matrix g({{2, 4, 6}}); + // Matrix g({{2, 4, 6}}); Matrix h({{7}, {9}, {11}}); Matrix i({{6, 1, 1}, {4, -2, 5}, {2, 8, 7}}); Matrix _i_row({{6, 1, 1}, {2, 8, 7}}); Matrix _i_col({{6, 1}, {4, 5}, {2, 7}}); Matrix _i_swap({{4, -2, 5}, {6, 1, 1}, {2, 8, 7}}); Matrix j({{6, 1, 1}, {4, -2, 5}, {2, -8, 7}}); + Matrix k({{6, 1, 1}, {4, -2, 5}, {2, 8, 7}}); + k.swapRows(0, 1); + std::array l = {{11, 18, -20, 5}}; + Matrix m( + {{-2, 1, 0, 0}, {1, -2, 1, 0}, {0, 1, -2, 1}, {0, 0, 1, -2}}); + std::array n = {{4, 7, 2, 5}}; + auto o = m * n; + std::cout << BLUE << "[ TESTS BEGINNING ]\n\n" << RESET; test(d + e, _de_add, "matrix addition"); test(d - e, _de_sub, "matrix subtraction"); test(b * c, _bc, "matrix multiplication"); @@ -44,29 +56,65 @@ int main() test(-e, _e_neg, "unary minus, negation"); test(identity(), _id, "identity construction"); test(determinant(f), -306, "determinant"); - test(dotProduct(g, h), 116, "dot product"); + // test(dotProduct(g, h), 116, "dot product"); test(removeRow(i, 1), _i_row, "remove row"); test(removeCol(i, 1), _i_col, "remove col"); - - Matrix k({{6, 1, 1}, {4, -2, 5}, {2, 8, 7}}); - k.swapRows(0, 1); test(k, _i_swap, "swap row"); test(j.findLargestInCol(1, 0), 2u, "find largest element in column"); + test(pNorm(l, 1), 54, "vector one norm"); + test(infNorm(l), 20, "vector infinity norm"); + test(oneNorm(i), 13, "matrix one norm"); + test(infNorm(i), 17, "matrix infinity norm"); + test(pNorm(jacobiIteration(m, o) - n, 2) < 1e-10, true, "jacobi iteration"); + + std::cout << BLUE << "\n[ TESTS COMPLETE ]\n\n" << RESET; +} + +int main() +{ + // runTests(); + + // Matrix U({{3, 5, -6, 4}, {0, 4, -6, 9}, {0, 0, 3, 11}, {0, 0, + // 0, -9}}); std::array x{4, 6, -7, 9}; auto B = U * x; + // + // std::cout << " U\n" << U << std::endl; + // std::cout << " b\n" << B << std::endl; + // std::cout << " Real x\n" << x << std::endl; + // std::cout << " Calculated x\n"; + // std::cout << U.backSub(B) << std::endl; + // + // std::cout << " v" << l << std::endl; + // std::cout << "1 Norm: " << pNorm(l, 1) << std::endl; + // std::cout << "2 Norm: " << pNorm(l, 2) << std::endl; + + // Matrix A({{-2, 1, 0, 0}, {1, -2, 1, 0}, {0, 1, -2, 1, 0}, {0, + // 0, 1, -2}}); Matrix A({{6, -1}, {2, 3}}); Matrix A( + // [](unsigned int const& i, unsigned int const& j) { return 1.0 / (i + j + // + 1.0); }); + + // auto eigval1 = inversePowerIteration(A, 1000u); + // std::cout << "A\n" << A << std::endl; + // std::cout << "Smallest Eigenvalue\n" << eigval1<< std::endl; + + // auto eigval = powerIteration(A, 1000u); + // std::cout << "A\n" << A << std::endl; + // std::cout << "Largest Eigenvalue\n" << eigval << std::endl; + + // std::cout << fivePointStencil() << std::endl; + // std::cout << ninePointStencil() << std::endl; + + // auto answer = solveNinePointStencil(0.0, 1.0, sin); + // auto finalMat = arrayToMat(answer); + // std::cout << "Answer in Matrix Form\n" << finalMat << std::endl; - // Matrix x({{3, -7, -2, 2}, {-3, 5, 1, 0}, {6, -4, 0, -5}, {-9, 5, -5, 12}}); - // //Matrix x({{2, 1, 1, 0}, {4, 3, 3, 1}, {8, 7, 9, 5}, {6, 7, 9, 8}}); - // auto [L, U, P] = x.luFactorize(); - // std::cout << "A\n" << x << std::endl; - // std::cout << "L\n" << L << std::endl << "U\n" << U << std::endl; - // std::cout << "PA\n" << P*x << std::endl; - // std::cout << "LU\n" << L*U << std::endl; + // Matrix A(1, 10); // random 5x5 with values from 0-10 + Matrix A({{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}); + auto [L, U, P] = A.luFactorize(); - Matrix l( - {{-2, 1, 0, 0}, {1, -2, 1, 0}, {0, 1, -2, 1}, {0, 0, 1, -2}}); - //auto foo = l.triDiagThomas({0, 1, 1, 1}, {-2, -2, -2, -2}, {1, 1, 1, 0}, {0.04, 0.04, 0.04, 0.04}); - auto foo = l.triDiagThomas({0.04, 0.04, 0.04, 0.04}); - for(auto && e:foo) - { - std::cout << e << " " << std::endl; - } + std::cout << " L\n" << L << std::endl; + std::cout << " U\n" << U << std::endl; + std::cout << " P\n" << P << std::endl << std::endl; + std::cout << " LU\n" << L * U << std::endl; + std::cout << " PA\n" << P * A << std::endl; } diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/manual_back_sub.md b/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/manual_back_sub.md new file mode 100644 index 0000000..defb649 --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/manual_back_sub.md @@ -0,0 +1,88 @@ +--- +title: Back Substitution +math: true +layout: default +--- + +{% include mathjax.html %} + + Table of Contents +# Back Substitution + +**Routine Name:** backSub + +**Author:** Philip Nelson + +**Language:** C++ + +## Description + +`backSub` solves a linear system \\(Ux=b\\) for \\(x\\) by back substitution where \\(U\\) is an upper triangular matrix. + +## Input + +`backSub(std::array b)` is called by a `Matrix of type `T` and size `MxM` and requires: + +* `std::array b` - a column vector `b` of type `T` and size `M` + +## Output + +`backSub` returns a `std::array` with the solution vector \\(x\\) + +## Code +{% highlight c++ %} +std::array backSub(std::array b) +{ + std::array x; + for (auto i = (int)M - 1; i >= 0; --i) + { + T sum = 0.0; + for (auto j = (unsigned int)i + 1; j < M; ++j) + { + sum += m[i][j] * x[j]; + } + x[i] = (b[i] - sum) / m[i][i]; + } + return x; +} +{% endhighlight %} + +## Example +{% highlight c++ %} +int main() +{ + Matrix U({ { 3, 5, -6, 4}, + { 0, 4, -6, 9}, + { 0, 0, 3, 11}, + { 0, 0, 0, -9} }); + + std::array x{4, 6, -7, 9}; + auto b = U * x; + + std::cout << " U\n" << U << std::endl; + std::cout << " b\n" << b << std::endl; + std::cout << " Real x\n" << x << std::endl; + std::cout << " Calculated x\n"; + std::cout << U.backSub(b) << std::endl; +} +{% endhighlight %} + +## Result +``` + U +| 3 5 -6 4 | +| 0 4 -6 9 | +| 0 0 3 11 | +| 0 0 0 -9 | + + b +[ 120 147 78 -81 ] + + Real x +[ 4 6 -7 9 ] + + Calculated x +[ 4 6 -7 9 ] +``` + +**Last Modification date:** 07 February 2018 diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/manual_forward_sub.md b/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/manual_forward_sub.md new file mode 100644 index 0000000..d224581 --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/manual_forward_sub.md @@ -0,0 +1,90 @@ +--- +title: Forward Substitution +math: true +layout: default +--- + +{% include mathjax.html %} + + Table of Contents +# Forward Substitution + +**Routine Name:** forwardSub + +**Author:** Philip Nelson + +**Language:** C++ + +## Description + +`forwardSub` solves a linear system \\(Ly=b\\) for \\(y\\) by forward substitution where \\(L\\) is a lower triangular matrix. + +## Input + +`forwardSub(std::array b)` is called by a `Matrix of type `T` and size `MxM` and requires: + +* `std::array b` - a column vector `b` of type `T` and size `M` + +## Output + +`forwardSub` returns a `std::array` with the solution vector \\(y\\) + +## Code +{% highlight c++ %} +std::array forwardSub(std::array b) +{ + std::array y; + for (auto i = 0u; i < N; ++i) + { + T sum = 0.0; + for (auto j = 0u; j < i; ++j) + { + sum += m[i][j] * y[j]; + } + y[i] = b[i] - sum; + } + return y; +} +{% endhighlight %} + +## Example +{% highlight c++ %} +int main() +{ + Matrix L({ { 1, 0, 0, 0}, + { 5, 1, 0, 0}, + { 4, -6, 1, 0}, + {-4, 5, -9, 1} }); + + std::array y{3, 5, -6, 8}; + auto b = L * y; + + std::cout << " L\n" << L << std::endl; + std::cout << " b\n" << b << std::endl; + std::cout << " Real y\n" << y << std::endl; + std::cout << " Calculated y\n"; + std::cout << L.forwardSub(b) << std::endl; +} +{% endhighlight %} + +## Result +``` + L +| 1 0 0 0 | +| 5 1 0 0 | +| 4 -6 1 0 | +| -4 5 -9 1 | + + b +[ 3 20 -24 75 ] + + Real y +[ 3 5 -6 8 ] + + Calculated y +[ 3 5 -6 8 ] + + +``` + +**Last Modification date:** 07 February 2018 diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/manual_gen_five_point_stencil.md b/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/manual_gen_five_point_stencil.md new file mode 100644 index 0000000..afb244e --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/manual_gen_five_point_stencil.md @@ -0,0 +1,68 @@ +--- +title: Five Point Stencil +math: true +layout: default +--- + +{% include mathjax.html %} + + Table of Contents +# Five Point Stencil + +**Routine Name:** fivePointStencil + +**Author:** Philip Nelson + +**Language:** C++ + +## Description + +`fivePointStencil` returns an \\(N^2\\)x\\(N^2\\) matrix for the five point stencil for the 2D finite difference method + +## Input + +`Matrix fivePointStencil()` takes no input but requires the size and template parameters \\(N\\) and \\(T\\) + +## Output + +A \\(N^2\\)x\\(N^2\\) matrix with the stencil + +## Code +{% highlight c++ %} +template +Matrix fivePointStencil() +{ + return Matrix([](int i, int j) { + if (i == j) return -4; + if (i + 1 == j) return i % 3 != 2 ? 1 : 0; + if (i == j + 1) return i % 3 != 0 ? 1 : 0; + if (i + 3 == j) return 1; + if (i == j + 3) return 1; + return 0; + + }); +} +{% endhighlight %} + +## Example +{% highlight c++ %} +int main() +{ + std::cout << fivePointStencil() << std::endl; +} +{% endhighlight %} + +## Result +``` +| -4 1 0 1 0 0 0 0 0 | +| 1 -4 1 0 1 0 0 0 0 | +| 0 1 -4 0 0 1 0 0 0 | +| 1 0 0 -4 1 0 1 0 0 | +| 0 1 0 1 -4 1 0 1 0 | +| 0 0 1 0 1 -4 0 0 1 | +| 0 0 0 1 0 0 -4 1 0 | +| 0 0 0 0 1 0 1 -4 1 | +| 0 0 0 0 0 1 0 1 -4 | +``` + +**Last Modification date:** 27 February 2018 diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/manual_gen_mesh.md b/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/manual_gen_mesh.md new file mode 100644 index 0000000..826a940 --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/manual_gen_mesh.md @@ -0,0 +1,64 @@ +--- +title: generateMesh +math: true +layout: default +--- + +{% include mathjax.html %} + + Table of Contents +# Generate Mesh + +**Routine Name:** generateMesh + +**Author:** Philip Nelson + +**Language:** C++ + +## Description + +`generateMesh(T a, T b)` generates a 2D mesh from `a` to `b` with `N` steps + +## Input + +`generateMesh(T a, T b)` requires: + +* `T a` - the start of the mesh +* `T b` - the end of the mesh +* `N` - the size + +## Output + +The mesh of \\(i\cdot j\\) + +## Code +{% highlight c++ %} +template +Matrix generateMesh(T a, T b) +{ + auto h = (b - a) / (N - 1); + + return Matrix( + [&](int i, int j) { return (a + (i + 1) * h) * (a + (j + 1) * h); }); +} +{% endhighlight %} + +## Example +{% highlight c++ %} +int main() +{ + std::cout << "mesh\n" << generateMesh(0, 1) << std::endl; +} +{% endhighlight %} + +## Result +``` +mesh +| 0.0278 0.0556 0.0833 0.111 0.139 | +| 0.0556 0.111 0.167 0.222 0.278 | +| 0.0833 0.167 0.25 0.333 0.417 | +| 0.111 0.222 0.333 0.444 0.556 | +| 0.139 0.278 0.417 0.556 0.694 | +``` + +**Last Modification date:** 27 February 2018 diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/manual_gen_nine_point_stencil.md b/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/manual_gen_nine_point_stencil.md new file mode 100644 index 0000000..47e13f8 --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/manual_gen_nine_point_stencil.md @@ -0,0 +1,72 @@ +--- +title: Nine Point Stencil +math: true +layout: default +--- + +{% include mathjax.html %} + + Table of Contents +# Nine Point Stencil + +**Routine Name:** ninePointStencil + +**Author:** Philip Nelson + +**Language:** C++ + +## Description + +`ninePointStencil` returns an \\(N^2\\)x\\(N^2\\) matrix for the nine point stencil for 2D finite difference method + +## Input + +`Matrix ninePointStencil()` takes no input but requires the size and template parameters \\(N\\) and \\(T\\) + +## Output + +A \\(N^2\\)x\\(N^2\\) matrix with the stencil + +## Code +{% highlight c++ %} +template +Matrix ninePointStencil() +{ + return Matrix([](int i, int j) { + if (i == j) return -20; + if (i + 1 == j) return i % 3 != 2 ? 4 : 0; + if (i == j + 1) return i % 3 != 0 ? 4 : 0; + if (i + 2 == j) return i % 3 != 0 ? 1 : 0; + if (i == j + 2) return i % 3 != 2 ? 1 : 0; + ; + if (i + 3 == j) return 4; + if (i == j + 3) return 4; + if (i + 4 == j) return i % 3 != 2 ? 1 : 0; + if (i == j + 4) return i % 3 != 0 ? 1 : 0; + ; + }); +} +{% endhighlight %} + +## Example +{% highlight c++ %} +int main() +{ + std::cout << ninePointStencil() << std::endl; +} +{% endhighlight %} + +## Result +``` +| -20 4 0 4 1 0 0 0 0 | +| 4 -20 4 1 4 1 0 0 0 | +| 0 4 -20 0 1 4 0 0 0 | +| 4 1 0 -20 4 0 4 1 0 | +| 1 4 1 4 -20 4 1 4 1 | +| 0 1 4 0 4 -20 0 1 4 | +| 0 0 0 4 1 0 -20 4 0 | +| 0 0 0 1 4 1 4 -20 4 | +| 0 0 0 0 1 4 0 4 -20 | +``` + +**Last Modification date:** 27 February 2018 diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/manual_init_b.md b/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/manual_init_b.md new file mode 100644 index 0000000..e804d0b --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/manual_init_b.md @@ -0,0 +1,66 @@ +--- +title: Initialize B for Mesh +math: true +layout: default +--- + +{% include mathjax.html %} + + Table of Contents +# Initialize B for Mesh + +**Routine Name:** initMeshB + +**Author:** Philip Nelson + +**Language:** C++ + +## Description + +`initMeshB` initializes the `b` matrix used for the 2D finite difference method + +## Input + +`initMeshB(Matrix& mesh, F f)` requires: + +* `Matrix` - the `N`x`N` mesh +* `F f` - the forcing function + +## Output + +The `b` matrix used to solve the 2D finite difference method + +## Code +{% highlight c++ %} +template +std::array initMeshB(Matrix& mesh, F f) +{ + std::array b; + for (auto i = 0u; i < N; ++i) + { + for (auto j = 0u; j < N; ++j) + { + b[i * N + j] = f(mesh[i][j]); + } + } + return b; +} +{% endhighlight %} + +## Example +{% highlight c++ %} +int main() +{ + auto mesh = generateMesh(0, 1); + auto b = initMeshB(mesh, sin); + std::cout << "b\n" << b << std::endl << std::endl; +} +{% endhighlight %} + +## Result +``` +b +[ 0.0278 0.0555 0.0832 0.111 0.138 0.0555 0.111 0.166 0.22 0.274 0.0832 0.166 0.247 0.327 0.405 0.111 0.22 0.327 0.43 0.527 0.138 0.274 0.405 0.527 0.64 ] +``` + +**Last Modification date:** <++>10 February 2018 diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/manual_inverse_power_iteration.md b/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/manual_inverse_power_iteration.md new file mode 100644 index 0000000..26b841d --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/manual_inverse_power_iteration.md @@ -0,0 +1,80 @@ +--- +title: Inverse Power Iteration +math: true +layout: default +--- + +{% include mathjax.html %} + + Table of Contents +# Power Iteration + +**Routine Name:** inversePowerIteration + +**Author:** Philip Nelson + +**Language:** C++ + +## Description + +`inversePowerIteration` calculates the smallest Eigenvalue of a `N`x`N` matrix. + +## Input + +`inversePowerIteration(Matrix const& A, unsigned int const& MAX)` requires: + +* `Matrix const& A` - an `N`x`N` matrix +* `unsigned int const& MAX` - the maximum number of iterations + +## Output + +The smallest Eigenvalue + +## Code +{% highlight c++ %} +template +T inversePowerIteration(Matrix & A, unsigned int const& MAX) +{ + std::array v; + for (auto&& e : v) + e = randDouble(0.0, 10.0); + + T lamda = 0; + for (auto i = 0u; i < MAX; ++i) + { + auto w = A.solveLinearSystemLU(v); + v = w / pNorm(w,2); + lamda = v*(A*v); + } + return lamda; +} +{% endhighlight %} + +## Example +{% highlight c++ %} +int main() +{ + Matrix A( + [](unsigned int const& i, unsigned int const& j) { return 1.0 / (i + j + 1.0); }); + + auto eigval = inversePowerIteration(A, 1000u); + std::cout << "A\n" << A << std::endl; + std::cout << "Smallest Eigenvalue\n" << eigval << std::endl; +} +{% endhighlight %} + +## Result +``` +A +| 1 0.5 0.333 0.25 0.2 | +| 0.5 0.333 0.25 0.2 0.167 | +| 0.333 0.25 0.2 0.167 0.143 | +| 0.25 0.2 0.167 0.143 0.125 | +| 0.2 0.167 0.143 0.125 0.111 | + +Smallest Eigenvalue +3.29e-06 + +``` + +**Last Modification date:** 27 February 2018 diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/manual_jacobi_iteration.md b/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/manual_jacobi_iteration.md new file mode 100644 index 0000000..1283b43 --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/manual_jacobi_iteration.md @@ -0,0 +1,108 @@ +--- +title: Jacobi Iteration +math: true +layout: default +--- + +{% include mathjax.html %} + + Table of Contents +# Jacobi Iteration + +**Routine Name:** jacobiIteration + +**Author:** Philip Nelson + +**Language:** C++ + +## Description + +`jacobiIteration` solves a linear system of equations by Jacobi Iteration. + +## Input + +`jacobiIteration(std::array const& b, unsigned int const& MAX)` is called by a `Matrix of type `T` and size `M`x`M` and requires: + +* `std::array b` - a column vector `b` of type `T` and size `M` +* `unsigned int MAX` - the max iterations (optional) + +## Output + +`jacobiIteration` returns a `std::array` with the solution vector `x` + +## Code +{% highlight c++ %} +std::array jacobiIteration(std::array const& b, unsigned int const& MAX = 1000) +{ + std::array zeros; + zeros.fill(0); + + std::array x = zeros; + + for (auto n = 0u; n < MAX; ++n) + { + auto x_n = zeros; + + for (auto i = 0u; i < M; ++i) + { + T sum = 0; + for (auto j = 0u; j < N; ++j) + { + if (j == i) + continue; + sum += m[i][j] * x[j]; + } + x_n[i] = (b[i] - sum) / m[i][i]; + } + + if (allclose(x, x_n, maceps().maceps)) + break; + + x = x_n; + } + + return x; +} +{% endhighlight %} + +## Example +{% highlight c++ %} +int main() +{ + Matrix A({ + {-2, 1, 0, 0}, + {1, -2, 1, 0}, + {0, 1, -2, 1}, + {0, 0, 1, -2} + }); + std::array x = {4, 7, 2, 5}; + auto b = A*x; + std::cout << " A\n" << A << std::endl; + std::cout << " b\n" << b << std::endl; + std::cout << " x\n" << x << std::endl; + + std::cout << "Calculated x\n"; + std::cout << A.jacobiIteration(b) << std::endl; +} +{% endhighlight %} + +## Result +``` + A +| -2 1 0 0 | +| 1 -2 1 0 | +| 0 1 -2 1 | +| 0 0 1 -2 | + + b +[ -1 -8 8 -8 ] + + x +[ 4 7 2 5 ] + +Calculated x +[ 4 7 2 5 ] + +``` + +**Last Modification date:** 07 February 2018 diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/manual_linear_solve_lu.md b/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/manual_linear_solve_lu.md new file mode 100644 index 0000000..aa756ba --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/manual_linear_solve_lu.md @@ -0,0 +1,83 @@ +--- +title: Linear Solver LU +math: true +layout: default +--- + +{% include mathjax.html %} + + Table of Contents +# Linear Solver by LU Factorization + +**Routine Name:** solveLinearSystemLU + +**Author:** Philip Nelson + +**Language:** C++ + +## Description + +`solveLinearSystemLU` solves a linear system of equations \\(Ax=b\\) by LU Factorization. The method used is: + +\\[LU=PA\\] +\\[LUx = Pb\\] +\\[Ly = Pb\\] +\\[Ux=y\\] + +## Input + +`solveLinearSystemLU(std::array b)` is called by a `Matrix` of type `T` and size `M`x`M` and requires: + +* `std::array b` - a column vector `b` of type `T` and size `M` + +## Output + +`solveLinearSystemLU` returns a `std::array` with the solution vector `x` + +## Code +{% highlight c++ %} +std::array solveLinearSystemLU(std::array b) +{ + auto[L, U, P] = luFactorize(); + auto y = L.forwardSub(P * b); + auto x = U.backSub(y); + return x; +} +{% endhighlight %} + +## Example +{% highlight c++ %} +int main() +{ + Matrix A(1, 10); // random 4x4 with values from 0-10 + std::array x = {4, 7, 2, 5}; + auto b = A*x; + std::cout << " A\n" << A << std::endl; + std::cout << " b\n" << b << std::endl; + std::cout << " x\n" << x << std::endl; + + std::cout << "Calculated x\n"; + std::cout << A.solveLinearSystemLU(b) << std::endl; +} +{% endhighlight %} + +## Result +``` + A +| 3 8 1 4 | +| 6 10 5 8 | +| 8 4 8 6 | +| 8 10 8 1 | + + b +[ 90 144 106 123 ] + + Real x +[ 4 7 2 5 ] + +Calculated x +[ 4 7 2 5 ] + +``` + +**Last Modification date:** 07 February 2018 diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/manual_lu_factorize.md b/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/manual_lu_factorize.md new file mode 100644 index 0000000..87e4174 --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/manual_lu_factorize.md @@ -0,0 +1,123 @@ +--- +title: LU Factorization +math: true +layout: default +--- + +{% include mathjax.html %} + + Table of Contents +# LU Factorization + +**Routine Name:** luFactorize + +**Author:** Philip Nelson + +**Language:** C++ + +## Description + +`luFactorize` turns a matrix into upper and lower triangular matrices such that \\(LU=PA\\) + +## Input + +`luFactorize()` takes not arguments but must be called by a `Matrix of type `T` and size `MxM` + +## Output + +`luFactorize` returns a `std::tuple, Matrix, Matrix>` with \\(L,\ U,\ P\\) matricies. Destructured returning is recommended. + +## Code +{% highlight c++ %} +std::tuple, Matrix, Matrix> luFactorize() +{ + auto I = identity(); + auto P = I; + Matrix L(0); + Matrix U(m); + for (auto j = 0u; j < N; ++j) // columns + { + auto largest = U.findLargestInCol(j, j); + if (largest != j) + { + L.swapRows(j, largest); + U.swapRows(j, largest); + P.swapRows(j, largest); + } + auto pivot = U[j][j]; + auto mod = I; + for (auto i = j + 1; i < N; ++i) // rows + { + mod[i][j] = -U[i][j] / pivot; + } + L = L + I - mod; + U = mod * U; + } + L = I + L; + return {L, U, P}; +} +{% endhighlight %} + +## Example +{% highlight c++ %} +int main() +{ + Matrix A(1, 10); // random 5x5 with values from 0-10 + auto [L, U, P] = A.luFactorize(); + + std::cout << " L\n" << L << std::endl; + std::cout << " U\n" << U << std::endl; + std::cout << " P\n" << P << std::endl << std::endl; + std::cout << " LU\n" << L*U << std::endl; + std::cout << " PA\n" << P*A << std::endl; +} +{% endhighlight %} + +## Result +``` + A +| 8 8 4 2 6 | +| 5 5 5 3 1 | +| 10 3 10 3 3 | +| 5 2 9 4 8 | +| 10 3 7 7 4 | + + L +| 1 0 0 0 0 | +| 0.8 1 0 0 0 | +| 0.5 0.08929 1 0 0 | +| 1 0 -0.6885 1 0 | +| 0.5 0.625 0.5738 0.05136 1 | + + U +| 10 3 10 3 3 | +| 0 5.6 -4 -0.4 3.6 | +| 0 0 4.357 2.536 6.179 | +| 0 0 0 5.746 5.254 | +| 0 0 0 0 -6.565 | + + P +| 0 0 1 0 0 | +| 1 0 0 0 0 | +| 0 0 0 1 0 | +| 0 0 0 0 1 | +| 0 1 0 0 0 | + + + LU +| 10 3 10 3 3 | +| 8 8 4 2 6 | +| 5 2 9 4 8 | +| 10 3 7 7 4 | +| 5 5 5 3 1 | + + PA +| 10 3 10 3 3 | +| 8 8 4 2 6 | +| 5 2 9 4 8 | +| 10 3 7 7 4 | +| 5 5 5 3 1 | + +``` + +**Last Modification date:** 07 February 2018 diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/manual_matrix_infinity_norm.md b/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/manual_matrix_infinity_norm.md new file mode 100644 index 0000000..35bd3f0 --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/manual_matrix_infinity_norm.md @@ -0,0 +1,71 @@ +--- +title: Matrix Infinity Norm +math: true +layout: default +--- + +{% include mathjax.html %} + + Table of Contents +# Matrix Infinity Norm + +**Routine Name:** infNorm + +**Author:** Philip Nelson + +**Language:** C++ + +## Description + +`infNorm` calculates the infinity norm of a matrix. The infinity norm, \\(\|A\|_{\infty }\\) is the largest absolute sum of the rows of \\(A\\). + +## Input + +`infNorm(Matrix& m)` requires: + +* `Matrix m` - an `M`x`N` matrix of type `T` + +## Output + +A double with the value of the infinity norm + +## Code +{% highlight c++ %} +template +T infNorm(Matrix& m) +{ + std::array rowSum; + for (auto i = 0u; i < N; ++i) + rowSum[i] = std::accumulate( + m.begin(i), m.end(i), 0, [](T sum, T elem) { return sum + std::abs(elem); }); + + return *std::max_element(rowSum.begin(), rowSum.end()); +} +{% endhighlight %} + +## Example +{% highlight c++ %} +int main() +{ + Matrix A({ + {6, 1, 1}, + {4, -2, 5}, + {2, 8, 7} + }); + + std::cout << " A\n" << A << std::endl; + std::cout << "Infinity Norm: " << infNorm(A) << std::endl; +} +{% endhighlight %} + +## Result +``` + A +| 6 1 1 | +| 4 -2 5 | +| 2 8 7 | + +Infinity Norm: 17 +``` + +**Last Modification date:** 10 February 2018 diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/manual_matrix_one_norm.md b/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/manual_matrix_one_norm.md new file mode 100644 index 0000000..b58cdaf --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/manual_matrix_one_norm.md @@ -0,0 +1,77 @@ +--- +title: Matrix One Norm +math: true +layout: default +--- + +{% include mathjax.html %} + + Table of Contents +# Matrix One Norm + +**Routine Name:** oneNorm + +**Author:** Philip Nelson + +**Language:** C++ + +## Description + +`oneNorm` calculates the one norm of a matrix. Give a matrix \\(A\\), the one norm, \\(\|A\|_{1}\\), is the maximum of the absolute column sums. + +## Input + + +`oneNorm(Matrix& m)` requires: + +* `Matrix m` - an `M`x`N` matrix of type `T` + +## Output + +A double with the desired one norm of the matrix. + +## Code +{% highlight c++ %} +template +T oneNorm(Matrix& m) +{ + std::array colSum; + for (auto j = 0u; j < M; ++j) + { + colSum[j] = 0; + for (auto i = 0u; i < N; ++i) + { + colSum[j] += std::abs(m[i][j]); + } + } + + return *std::max_element(colSum.begin(), colSum.end()); +} +{% endhighlight %} + +## Example +{% highlight c++ %} +int main() +{ +Matrix A({ + {6, 1, 1}, + {4, -2, 5}, + {2, 8, 7} + }); + +std::cout << " A\n" << A << std::endl; +std::cout << "One Norm: " << oneNorm(A) << std::endl; +} +{% endhighlight %} + +## Result +``` + A +| 6 1 1 | +| 4 -2 5 | +| 2 8 7 | + +One Norm: 13 +``` + +**Last Modification date:** 10 February 2018 diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/manual_power_iteration.md b/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/manual_power_iteration.md new file mode 100644 index 0000000..64c7c01 --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/manual_power_iteration.md @@ -0,0 +1,81 @@ +--- +title: Power Iteration +math: true +layout: default +--- + +{% include mathjax.html %} + + Table of Contents +# Power Iteration + +**Routine Name:** powerIteration + +**Author:** Philip Nelson + +**Language:** C++ + +## Description + +`powerIteration` calculates the largest Eigenvalue of a `N`x`N` matrix. + +## Input + + +`powerIteration(Matrix const& A, unsigned int const& MAX)` requires: + +* `Matrix const& A` - an `N`x`N` matrix +* `unsigned int const& MAX` - the maximum number of iterations + +## Output + +The largest Eigenvalue + +## Code +{% highlight c++ %} +template +T powerIteration(Matrix const& A, unsigned int const& MAX) +{ + std::array b_k; + + for (auto&& e : b_k) + e = randDouble(0.0, 10.0); + + for (auto i = 0u; i < MAX; ++i) + { + auto Ab_k = A * b_k; + auto norm = pNorm(Ab_k, 2); + b_k = Ab_k / norm; + } + return pNorm(A * b_k, 2); +} +{% endhighlight %} + +## Example +{% highlight c++ %} +int main() +{ + Matrix A( + [](unsigned int const& i, unsigned int const& j) { return 1.0 / (i + j + 1.0); }); + + auto eigval = powerIteration(A, 1000u); + std::cout << "A\n" << A << std::endl; + std::cout << "Largest Eigenvalue\n" << eigval << std::endl; +} +{% endhighlight %} + +## Result +``` +A +| 1 0.5 0.333 0.25 0.2 | +| 0.5 0.333 0.25 0.2 0.167 | +| 0.333 0.25 0.2 0.167 0.143 | +| 0.25 0.2 0.167 0.143 0.125 | +| 0.2 0.167 0.143 0.125 0.111 | + +Largest Eigenvalue +1.57 + +``` + +**Last Modification date:** 27 February 2018 diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/manual_solve_five_point_stencil.md b/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/manual_solve_five_point_stencil.md new file mode 100644 index 0000000..cccee13 --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/manual_solve_five_point_stencil.md @@ -0,0 +1,98 @@ +--- +title: 5-Point Stencil +math: true +layout: default +--- + +{% include mathjax.html %} + + Table of Contents +# Five Point Stencil + +**Routine Name:** solveFivePointStencil + +**Author:** Philip Nelson + +**Language:** C++ + +## Description + +`fivePointStencil` solves Poisson's Equation for an arbitrary mesh size and forcing function + +## Input + +`solveFivePointStencil(T a, T b, F f)` requires: + +* `T a` - the start boundary +* `T b` - the end boundary +* `F f` - the forcing function +* `T` - template param for the type +* `N` - template param for the dimension of the mesh + +## Output + +An array with the approximated solution + +## Code +{% highlight c++ %} +template +auto solveFivePointStencil(T a, T b, F f) +{ + auto mesh = generateMesh(a, b); + auto bv = initMeshB(mesh, f); + auto stencil = fivePointStencil(); + auto res = stencil.solveLinearSystemLU(bv); + return res; +} +{% endhighlight %} + +solveFivePointStencil relies on: +[fivePointStencil](./manual_gen_five_point_stencil)| +[generateMesh](./manual_gen_mesh)| +[initMeshB](./manual_init_b)| +[solveLinearSystemLU](./manual_linear_solve_lu)| + +## Example +{% highlight c++ %} +int main() +{ + auto answer = solveFivePointStencil(0.0, 1.0, sin); + auto finalMat = arrayToMat(answer); + std::cout << "Answer in Matrix Form\n" << finalMat << std::endl; +} +{% endhighlight %} + +## Result +``` +mesh +| 0.0625 0.125 0.188 | +| 0.125 0.25 0.375 | +| 0.188 0.375 0.562 | + + +stencil +| -4 1 0 1 0 0 0 0 0 | +| 1 -4 1 0 1 0 0 0 0 | +| 0 1 -4 0 0 1 0 0 0 | +| 1 0 0 -4 1 0 1 0 0 | +| 0 1 0 1 -4 1 0 1 0 | +| 0 0 1 0 1 -4 0 0 1 | +| 0 0 0 1 0 0 -4 1 0 | +| 0 0 0 0 1 0 1 -4 1 | +| 0 0 0 0 0 1 0 1 -4 | + + +bv +[ 0.0625 0.125 0.186 0.125 0.247 0.366 0.186 0.366 0.533 ] + + +result +[ -0.097 -0.163 -0.154 -0.163 -0.276 -0.266 -0.154 -0.266 -0.266 ] + +Answer in Matrix Form +| -0.097 -0.163 -0.154 | +| -0.163 -0.276 -0.266 | +| -0.154 -0.266 -0.266 | +``` + +**Last Modification date:** 28 February 2018 diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/manual_solve_nine_point_stencil.md b/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/manual_solve_nine_point_stencil.md new file mode 100644 index 0000000..ae4dea1 --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/manual_solve_nine_point_stencil.md @@ -0,0 +1,97 @@ +--- +title: 9-Point Stencil +math: true +layout: default +--- + +{% include mathjax.html %} + + Table of Contents +# Nine Point Stencil + +**Routine Name:** solveNinePointStencil + +**Author:** Philip Nelson + +**Language:** C++ + +## Description + +`solveNinePointStencil` solves Poisson's Equation for an arbitrary mesh size and forcing function + +## Input + +`solveNinePointStencil(T a, T b, F f)` requires: + +* `T a` - the start boundary +* `T b` - the end boundary +* `F f` - the forcing function +* `T` - template param for the type +* `N` - template param for the dimension of the mesh + +## Output + +An array with the approximated solution + +## Code +{% highlight c++ %} +template +auto solveFivePointStencil(T a, T b, F f) +{ + auto mesh = generateMesh(a, b); + auto bv = initMeshB(mesh, f); + auto stencil = ninePointStencil(); + auto res = stencil.solveLinearSystemLU(bv); + return res; +} +{% endhighlight %} + +solveFivePointStencil relies on: +[ninePointStencil](./manual_gen_nine_point_stencil)| +[generateMesh](./manual_gen_mesh)| +[initMeshB](./manual_init_b)| +[solveLinearSystemLU](./manual_linear_solve_lu)| + +## Example +{% highlight c++ %} +int main() +{ + auto answer = solveNinePointStencil(0.0, 1.0, sin); + auto finalMat = arrayToMat(answer); + std::cout << "Answer in Matrix Form\n" << finalMat << std::endl; +} +{% endhighlight %} + +## Result +``` +| 0.0625 0.125 0.188 | +| 0.125 0.25 0.375 | +| 0.188 0.375 0.562 | + + +stencil +| -20 4 0 4 1 0 0 0 0 | +| 4 -20 4 1 4 1 0 0 0 | +| 0 4 -20 0 1 4 0 0 0 | +| 4 1 0 -20 4 0 4 1 0 | +| 1 4 1 4 -20 4 1 4 1 | +| 0 1 4 0 4 -20 0 1 4 | +| 0 0 0 4 1 0 -20 4 0 | +| 0 0 0 1 4 1 4 -20 4 | +| 0 0 0 0 1 4 0 4 -20 | + + +bv +[ 0.0625 0.125 0.186 0.125 0.247 0.366 0.186 0.366 0.533 ] + + +result +[ -0.0169 -0.0284 -0.0267 -0.0284 -0.0483 -0.0466 -0.0267 -0.0466 -0.0477 ] + +Answer in Matrix Form +| -0.0169 -0.0284 -0.0267 | +| -0.0284 -0.0483 -0.0466 | +| -0.0267 -0.0466 -0.0477 | +``` + +**Last Modification date:** 28 February 2018 diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/manual_thomas_algorithm.md b/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/manual_thomas_algorithm.md new file mode 100644 index 0000000..effa19a --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/manual_thomas_algorithm.md @@ -0,0 +1,101 @@ +--- +title: Thomas Algorithm +math: true +layout: default +--- + +{% include mathjax.html %} + + Table of Contents +# Thomas Algorithm + +**Routine Name:** triDiagThomas + +**Author:** Philip Nelson + +**Language:** C++ + +## Description + +`triDiagThomas` solves a tridiagonal linear system of equations using the Thomas Algorithm. + +## Input + +`triDiagThomas(std::array const& d)` is called by a `Matrix` of type `T` and size `M`x`M` and requires: + + +* `std::array d` - an array of type `T` and size `M` (d vector) + +## Output + +`triDiagThomas` returns a `std::array` with the solution vector `x` + +## Code +{% highlight c++ %} +std::array triDiagThomas(std::array const& a, + std::array const& b, + std::array const& c, + std::array const& d) +{ + std::array cp, dp, x; + cp[0] = c[0] / b[0]; + dp[0] = d[0] / b[0]; + for (auto i = 1u; i < N; ++i) + { + double bottom = (b[i] - (a[i] * cp[i - 1])); + cp[i] = c[i] / bottom; + dp[i] = (d[i] - (a[i] * dp[i - 1])) / bottom; + } + + x[N - 1] = dp[N - 1]; + + for (auto i = (int)N - 2; i >= 0; --i) + { + x[i] = dp[i] - cp[i] * x[i + 1]; + } + return x; +} +{% endhighlight %} + +## Example +{% highlight c++ %} +int main() +{ + Matrix A({ + {-2, 1, 0, 0}, + { 1, -2, 1, 0}, + { 0, 1, -2, 1}, + { 0, 0, 1, -2} + }); + std::array x = {4, -3, 5, 1}; + + auto d = A * x; + auto testx = A.triDiagThomas(d); + + std::cout << " A\n" << A << std::endl; + std::cout << " d\n" << d << std::endl; + std::cout << " Real x\n" << x << std::endl; + std::cout << " Calculated x\n" << testx << std::endl; +} +{% endhighlight %} + +## Result +``` + A +| -2 1 0 0 | +| 1 -2 1 0 | +| 0 1 -2 1 | +| 0 0 1 -2 | + + d +[ -11 15 -12 3 ] + + Real x +[ 4 -3 5 1 ] + + Calculated x +[ 4 -3 5 1 ] + +``` + +**Last Modification date:** 07 February 2018 diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/manual_vector_infinity_norm.md b/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/manual_vector_infinity_norm.md new file mode 100644 index 0000000..1ce94d3 --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/manual_vector_infinity_norm.md @@ -0,0 +1,63 @@ +--- +title: Infinity Norm +math: true +layout: default +--- + +{% include mathjax.html %} + + Table of Contents +# Infinity Norm + +**Routine Name:** infNorm + +**Author:** Philip Nelson + +**Language:** C++ + +## Description + +`infNorm` calculates the infinity norm of a vector which is also the maximum element + +## Input + +`infNorm(std::array v)` requires: + +* `std::array v` - a column vector of type `T` and size `M` + +## Output + +A double with the value of the infinity norm + +## Code +{% highlight c++ %} +template +double infNorm(std::array v) +{ + T max = std::abs(v[0]); + for (auto&& x : v) + max = std::max(max, std::abs(x)); + return max; +} +{% endhighlight %} + +## Example +{% highlight c++ %} +int main() +{ + std::array v = {11, 18, -20, 5}; + + std::cout << " v\n" << v << std::endl; + std::cout << "Infinity Norm: " << infNorm(v) << std::endl; +} +{% endhighlight %} + +## Result +``` + v +[ 11 18 -20 5 ] + +Infinity Norm: 20 +``` + +**Last Modification date:** 07 February 2018 diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/manual_vector_pnorm.md b/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/manual_vector_pnorm.md new file mode 100644 index 0000000..aa93cd3 --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/manual_vector_pnorm.md @@ -0,0 +1,66 @@ +--- +title: P Norm +math: true +layout: default +--- + +{% include mathjax.html %} + + Table of Contents +# P Norm + +**Routine Name:** pNorm + +**Author:** Philip Nelson + +**Language:** C++ + +## Description + +`pNorm` calculates the p norm of a vector + +## Input + +`pNorm(std::array v, unsigned int const& p)` requires: + +* `std::array v` - a column vector of type `T` and size `M` +* `unsigned int p` - the p norm that is desired + +## Output + +A double with the desired norm value. + +## Code +{% highlight c++ %} +template +double pNorm(std::array v, unsigned int const& p) +{ + double sum = 0.0; + for (auto&& x : v) + sum += std::pow(std::abs(x), p); + return std::pow(sum, 1.0 / p); +} +{% endhighlight %} + +## Example +{% highlight c++ %} +int main() +{ + std::array v = {11, 18, -20, 5}; + + std::cout << " v\n" << v << std::endl; + std::cout << "1 Norm: " << pNorm(v, 1) << std::endl; + std::cout << "2 Norm: " << pNorm(v, 2) << std::endl; +} +{% endhighlight %} + +## Result +``` + v +[ 11 18 -20 5 ] + +1 Norm: 54 +2 Norm: 29.5 +``` + +**Last Modification date:** 07 February 2018 diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/matrix.hpp b/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/matrix.hpp index a36db8e..7a7526f 100644 --- a/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/matrix.hpp +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/matrix.hpp @@ -5,18 +5,17 @@ #include "random.hpp" #include #include +#include +#include #include #include template class Matrix; -// template -// using filler = std::function; - /* returns an NxN identity matrix */ template -Matrix identity() +static Matrix identity() { Matrix matrix(0); for (auto i = 0u; i < N; ++i) @@ -30,15 +29,24 @@ template class Matrix { public: + using filler = std::function; /* Default Creation */ Matrix() {} + /* Filled Creation */ + Matrix(filler f) + { + for (auto i = 0u; i < M; ++i) + for (auto j = 0u; j < N; ++j) + m[i][j] = f(i, j); + } + /* Random Creation */ Matrix(int start, int end) { for (auto i = 0u; i < M; ++i) for (auto j = 0u; j < N; ++j) - m[i][j] = rand(start, end); + m[i][j] = randInt(start, end); } /* Fill With n */ @@ -83,22 +91,26 @@ class Matrix T get(unsigned int const& i, unsigned int const& j) const { return m[i][j]; } - void set(unsigned int const& i, unsigned int const& j, T const& val) { m[i][j] = val; } + void set(unsigned int const& i, unsigned int const& j, T const& val) + { + m[i][j] = val; + } std::array& operator[](int x) { return m[x]; } + auto begin(unsigned int n) { return m[n].begin(); } + + auto end(unsigned int n) { return m[n].end(); } + /* Swap rows r1 and r2 */ void swapRows(unsigned int const& r1, unsigned int const& r2) { - for (auto i = 0u; i < N; ++i) - { - std::swap(m[r1][i], m[r2][i]); - } - // return this; + std::swap(m[r1], m[r2]); } /* return the absolute largest element of a col starting at a given row */ - unsigned int findLargestInCol(unsigned int const& col, unsigned int const& row = 0) + unsigned int findLargestInCol(unsigned int const& col, + unsigned int const& row = 0) { T max = row; for (auto i = row + 1; i < M; ++i) @@ -115,15 +127,13 @@ class Matrix std::swap(m[j][i], m[i][j]); } - /* calculate the lower and upper triangles */ + /* calculate the lower and upper triangles with a permutation matrix*/ std::tuple, Matrix, Matrix> luFactorize() { auto I = identity(); - auto P = identity(); - P.transpose(); + auto P = I; Matrix L(0); Matrix U(m); - std::vector> swaps; for (auto j = 0u; j < N; ++j) // columns { auto largest = U.findLargestInCol(j, j); @@ -132,42 +142,80 @@ class Matrix L.swapRows(j, largest); U.swapRows(j, largest); P.swapRows(j, largest); - swaps.push_back({j, largest}); } auto pivot = U[j][j]; - auto mod = identity(); + auto mod = I; for (auto i = j + 1; i < N; ++i) // rows { mod[i][j] = -U[i][j] / pivot; } - L = -(mod - I) + L; + L = L + I - mod; U = mod * U; } L = I + L; return {L, U, P}; } - std::array triDiagThomas(std::array const& a, - std::array const& b, - std::array const& c, - std::array const& d) + std::array backSub(std::array b) { - std::array c_s, d_s, f; - c_s[0] = c[0] / b[0]; - d_s[0] = d[0] / b[0]; - for (auto i = 1u; i < M; ++i) + std::array x; + for (auto i = (int)M - 1; i >= 0; --i) { - auto bmcsta = 1.0 / (b[i] - c_s[i - 1] * a[i]); - c_s[i] = c[i] * bmcsta; - d_s[i] = (d[i] - d_s[i - 1] * a[i]) * bmcsta; + T sum = 0.0; + for (auto j = (unsigned int)i + 1; j < M; ++j) + { + sum += m[i][j] * x[j]; + } + x[i] = (b[i] - sum) / m[i][i]; } + return x; + } - f[M - 1] = d_s[M - 1]; - for (auto i = M - 2; i-- > 0;) + std::array forwardSub(std::array b) + { + std::array y; + for (auto i = 0u; i < N; ++i) { - f[i] = d_s[i] - c_s[i] * d[i + 1]; + T sum = 0.0; + for (auto j = 0u; j < i; ++j) + { + sum += m[i][j] * y[j]; + } + y[i] = b[i] - sum; + } + return y; + } + + std::array solveLinearSystemLU(std::array b) + { + auto [L, U, P] = luFactorize(); + auto y = L.forwardSub(P * b); + auto x = U.backSub(y); + return x; + } + + static std::array triDiagThomas(std::array const& a, + std::array const& b, + std::array const& c, + std::array const& d) + { + std::array cp, dp, x; + cp[0] = c[0] / b[0]; + dp[0] = d[0] / b[0]; + for (auto i = 1u; i < N; ++i) + { + double bottom = (b[i] - (a[i] * cp[i - 1])); + cp[i] = c[i] / bottom; + dp[i] = (d[i] - (a[i] * dp[i - 1])) / bottom; } - return f; + + x[N - 1] = dp[N - 1]; + + for (auto i = (int)N - 2; i >= 0; --i) + { + x[i] = dp[i] - cp[i] * x[i + 1]; + } + return x; } std::array triDiagThomas(std::array const& d) @@ -186,19 +234,19 @@ class Matrix b[M - 1] = m[M - 1][M - 1]; c[M - 1] = 0; - for (auto e : a) - std::cout << e << " "; - std::cout << std::endl; - for (auto e : b) - std::cout << e << " "; - std::cout << std::endl; - for (auto e : c) - std::cout << e << " "; - std::cout << std::endl; - return triDiagThomas(a, b, c, d); } + Matrix blockToMatrix( + Matrix, N * N, N * N> block) + { + Matrix res([&](int i, int j) { + return block[i / N][j / N][i - (N * (i / N))][j - (N * (j / N))]; + }); + + return res; + } + private: std::array, M> m; }; diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/matrix_util.hpp b/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/matrix_util.hpp index 383859b..7829362 100644 --- a/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/matrix_util.hpp +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/matrix_util.hpp @@ -1,10 +1,23 @@ #ifndef MATRIX_UTIL_HPP #define MATRIX_UTIL_HPP +#include "../machineEpsilon/maceps.hpp" +#include "random.hpp" +#include "vector_util.hpp" +#include +#include #include #include +#include #include +struct Point +{ + Point(double x, double y) : x(x), y(y) {} + double x; + double y; +}; + template class Matrix; @@ -87,71 +100,117 @@ T determinant(Matrix const& a) return det; } -/* Addition */ +/* Addition and Subtraction */ +#define matrix_add_subtract(op) \ + template \ + Matrix operator op( \ + Matrix const& a, Matrix const& b) \ + { \ + R matrix[M][N]; \ + for (auto i = 0u; i < M; ++i) \ + { \ + for (auto j = 0u; j < N; ++j) \ + { \ + matrix[i][j] = a.get(i, j) op b.get(i, j); \ + } \ + } \ + return matrix; \ + } + +matrix_add_subtract(+) matrix_add_subtract(-) + + /* Multiplication Matrix * Matrix*/ + template + Matrix operator*(Matrix const& a, Matrix const& b) +{ + R matrix[M][O]; + for (auto i = 0u; i < M; ++i) + { + for (auto j = 0u; j < O; ++j) + { + matrix[i][j] = dotProduct(a, b, i, j); + } + } + return matrix; +} + +/* Multiplication Matrix * Array */ template -Matrix operator+(Matrix const& a, Matrix const& b) + typename R = decltype(T() * U())> +std::array operator*(Matrix const& a, std::array const& b) { - R matrix[M][N]; + std::array matrix; for (auto i = 0u; i < M; ++i) { + R sum = 0; for (auto j = 0u; j < N; ++j) { - matrix[i][j] = a.get(i, j) + b.get(i, j); + sum += a.get(i, j) * b[j]; } + matrix[i] = sum; } return matrix; } -/* Subtraction */ +/* Multiplication Array * Matrix */ template -Matrix operator-(Matrix const& a, Matrix const& b) + typename R = decltype(T() * U())> +std::array operator*(std::array const& b, Matrix const& a) { - R matrix[M][N]; + std::array matrix; for (auto i = 0u; i < M; ++i) { + R sum = 0; for (auto j = 0u; j < N; ++j) { - matrix[i][j] = a.get(i, j) - b.get(i, j); + sum += a.get(i, j) * b[j]; } + matrix[i] = sum; } return matrix; } -/* Multiplication */ +/* Scalar Multiplication a * scalar */ template -Matrix operator*(Matrix const& a, Matrix const& b) +Matrix operator*(Matrix const& a, U const& scalar) { - R matrix[M][O]; + R matrix[M][N]; for (auto i = 0u; i < M; ++i) { - for (auto j = 0u; j < O; ++j) + for (auto j = 0u; j < N; ++j) { - matrix[i][j] = dotProduct(a, b, i, j); + matrix[i][j] = (a.get(i, j) * scalar); } } return matrix; } -/* Scalar Multiplication a*scalar */ +/* Scalar Multiplication scalar*a */ template -Matrix operator*(Matrix const& a, U scalar) +Matrix operator*(U const& scalar, Matrix const& a) { R matrix[M][N]; for (auto i = 0u; i < M; ++i) @@ -164,20 +223,20 @@ Matrix operator*(Matrix const& a, U scalar) return matrix; } -/* Scalar Multiplication scalar*a */ +/* Scalar Division a/scalar */ template -Matrix operator*(U scalar, Matrix const& a) + typename R = decltype(T() / U())> +Matrix operator/(Matrix const& a, U const& scalar) { R matrix[M][N]; for (auto i = 0u; i < M; ++i) { for (auto j = 0u; j < N; ++j) { - matrix[i][j] = (a.get(i, j) * scalar); + matrix[i][j] = (a.get(i, j) / scalar); } } return matrix; @@ -225,11 +284,206 @@ std::ostream& operator<<(std::ostream& o, Matrix const& m) o << "| "; for (auto j = 0u; j < N; ++j) { - o << std::setw(9) << std::setprecision(4) << std::setfill(' ') << m.get(i, j); + o << std::setw(10) << std::setprecision(3) << std::setfill(' ') + << m.get(i, j); } o << " |" << std::endl; } return o; } +template +T oneNorm(Matrix& m) +{ + std::array colSum; + for (auto j = 0u; j < M; ++j) + { + colSum[j] = 0; + for (auto i = 0u; i < N; ++i) + { + colSum[j] += std::abs(m[i][j]); + } + } + + return *std::max_element(colSum.begin(), colSum.end()); +} + +template +T infNorm(Matrix& m) +{ + std::array rowSum; + for (auto i = 0u; i < N; ++i) + rowSum[i] = std::accumulate(m.begin(i), m.end(i), 0, [](T sum, T elem) { + return sum + std::abs(elem); + }); + + return *std::max_element(rowSum.begin(), rowSum.end()); +} + +template +T powerIteration(Matrix const& A, unsigned int const& MAX) +{ + std::array b_k; + + for (auto&& e : b_k) + e = randDouble(0.0, 10.0); + + for (auto i = 0u; i < MAX; ++i) + { + auto Ab_k = A * b_k; + auto norm = pNorm(Ab_k, 2); + b_k = Ab_k / norm; + } + return pNorm(A * b_k, 2); +} + +template +T inversePowerIteration(Matrix& A, unsigned int const& MAX) +{ + std::array v; + for (auto&& e : v) + e = randDouble(0.0, 10.0); + + T lamda = 0; + for (auto i = 0u; i < MAX; ++i) + { + auto w = A.solveLinearSystemLU(v); + v = w / pNorm(w, 2); + lamda = v * (A * v); + } + return lamda; +} + +template +Matrix fivePointStencil() +{ + return Matrix([](int i, int j) { + if (i == j) return -4; + if (i + 1 == j) return i % 3 != 2 ? 1 : 0; + if (i == j + 1) return i % 3 != 0 ? 1 : 0; + if (i + 3 == j) return 1; + if (i == j + 3) return 1; + return 0; + }); + // return { + // {-4, 1, 0, 1, 0, 0, 0, 0, 0}, + // {1, -4, 1, 0, 1, 0, 0, 0, 0}, + // {0, 1, -4, 0, 0, 1, 0, 0, 0}, + // {1, 0, 0, -4, 1, 0, 1, 0, 0}, + // {0, 1, 0, 1, -4, 1, 0, 1, 0}, + // {0, 0, 1, 0, 1, -4, 0, 0, 1}, + // {0, 0, 0, 1, 0, 0, -4, 1, 0}, + // {0, 0, 0, 0, 1, 0, 1, -4, 1}, + // {0, 0, 0, 0, 0, 1, 0, 1, -4} + // }; +} + +template +Matrix ninePointStencil() +{ + return Matrix([](int i, int j) { + if (i == j) return -20; + if (i + 1 == j) return i % 3 != 2 ? 4 : 0; + if (i == j + 1) return i % 3 != 0 ? 4 : 0; + if (i + 2 == j) return i % 3 != 0 ? 1 : 0; + if (i == j + 2) return i % 3 != 2 ? 1 : 0; + ; + if (i + 3 == j) return 4; + if (i == j + 3) return 4; + if (i + 4 == j) return i % 3 != 2 ? 1 : 0; + if (i == j + 4) return i % 3 != 0 ? 1 : 0; + ; + }); + // return { + // {-4, 1, 0, 1, 0, 0, 0, 0, 0}, + // {1, -4, 1, 0, 1, 0, 0, 0, 0}, + // {0, 1, -4, 0, 0, 1, 0, 0, 0}, + // {1, 0, 0, -4, 1, 0, 1, 0, 0}, + // {0, 1, 0, 1, -4, 1, 0, 1, 0}, + // {0, 0, 1, 0, 1, -4, 0, 0, 1}, + // {0, 0, 0, 1, 0, 0, -4, 1, 0}, + // {0, 0, 0, 0, 1, 0, 1, -4, 1}, + // {0, 0, 0, 0, 0, 1, 0, 1, -4} + // }; +} + +template +Matrix generateMesh(T a, T b) +{ + auto h = (b - a) / (N - 1); + + return Matrix( + [&](int i, int j) { return (a + (i + 1) * h) * (a + (j + 1) * h); }); +} + +template +std::array initMeshB(Matrix& mesh, F f) +{ + std::array b; + for (auto i = 0u; i < N; ++i) + { + for (auto j = 0u; j < N; ++j) + { + b[i * N + j] = f(mesh[i][j]); + } + } + return b; +} + +template +double conditionNumber(Matrix m) +{ + auto max = powerIteration(m, 1000u); + auto min = inversePowerIteration(m, 1000u); + + return max / min; +} + +template +Matrix arrayToMat(std::array a) +{ + return Matrix( + [&](int i, int j) { return a[i * int(std::sqrt(N)) + j]; }); +} + +template +auto solveFivePointStencil(T a, T b, F f) +{ + auto mesh = generateMesh(a, b); + auto bv = initMeshB(mesh, f); + auto stencil = fivePointStencil(); + auto res = stencil.solveLinearSystemLU(bv); + return res; +} + +template +auto solveNinePointStencil(T a, T b, F f) +{ + auto mesh = generateMesh(a, b); + auto bv = initMeshB(mesh, f); + auto stencil = ninePointStencil(); + auto res = stencil.solveLinearSystemLU(bv); + return res; +} + +template +Matrix makeNDiag(std::vector const& entries, + unsigned int const& offset) +{ + return Matrix( + [=](const unsigned int& row, const unsigned int& col) -> T { + const int start = row + offset; + const int end = start + entries.size() - 1; + + if ((int)col >= start && (int)col <= end) + { + return entries[(int)col - row - offset]; + } + else + { + return 0; + } + }); +} + #endif diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/random.hpp b/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/random.hpp index 290f968..74e6bfe 100644 --- a/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/random.hpp +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/random.hpp @@ -3,7 +3,7 @@ #include -int rand(int low, int high) +int randInt(int low, int high) { static std::random_device rd; static std::mt19937 mt(rd()); @@ -11,4 +11,12 @@ int rand(int low, int high) return dist(mt); } +double randDouble(double low, double high) +{ + static std::random_device rd; + static std::mt19937 mt(rd()); + std::uniform_real_distribution<> dist(low, high); + return dist(mt); +} + #endif diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/vector_util.hpp b/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/vector_util.hpp new file mode 100644 index 0000000..380ff4e --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/matrix/vector_util.hpp @@ -0,0 +1,141 @@ +#ifndef VECTOR_UTILL_HPP +#define VECTOR_UTILL_HPP + +#include "random.hpp" +#include +#include +#include +#include +#include + +template +std::array initRandom(int start, int end) +{ + std::array a; + for (auto i = 0u; i < N; ++i) + a[i] = randDouble(start, end); + + return a; +} + +template +bool allclose(std::array a, std::array b, double tol) +{ + for (auto i = 0u; i < M; ++i) + if (std::abs(a[i] - b[i]) > tol) return false; + return true; +} + +template +double pNorm(std::array v, unsigned int const& p) +{ + double sum = 0.0; + for (auto&& x : v) + sum += std::pow(std::abs(x), p); + return std::pow(sum, 1.0 / p); +} + +template +double infNorm(std::array v) +{ + T max = std::abs(v[0]); + for (auto&& x : v) + max = std::max(max, std::abs(x)); + return max; +} + +template +bool operator==(std::array a, std::array b) +{ + for (auto i = 0u; i < N; ++i) + if (a[i] != b[i]) return false; + return true; +} + +#define vector_add_subtract(op) \ + template \ + std::array operator op( \ + std::array const& a, std::array const& b) \ + { \ + std::array result; \ + for (auto i = 0u; i < N; ++i) \ + result[i] = a[i] op b[i]; \ + return result; \ + } + +vector_add_subtract(+) vector_add_subtract(-) + +#define vector_add_subtract_scalar(op) \ + template \ + R operator op(std::array const& a, U const& b) \ + { \ + R result = 0; \ + for (auto i = 0u; i < N; ++i) \ + result += a[i] op b; \ + return result; \ + } + + vector_add_subtract_scalar(+) vector_add_subtract_scalar(-) + +#define vector_multiply_divide_scalar(op) \ + template \ + std::array operator op(std::array const& a, U const& b) \ + { \ + std::array result; \ + for (auto i = 0u; i < N; ++i) \ + result[i] = a[i] op b; \ + return result; \ + } + + vector_multiply_divide_scalar(*) vector_multiply_divide_scalar(/) + +#define scalar_multiply_divide_vector(op) \ + template \ + std::array operator op(U const& b, std::array const& a) \ + { \ + std::array result; \ + for (auto i = 0u; i < N; ++i) \ + result[i] = b op a[i]; \ + return result; \ + } + + scalar_multiply_divide_vector(*) scalar_multiply_divide_vector(/) + + template + R operator*(std::array a, std::array b) +{ + R result = 0; + for (auto i = 0u; i < N; ++i) + result += a[i] * b[i]; + return result; +} + +// template +// std::ostream& operator<<(std::ostream& o, std::vector const& a) +//{ +// o << "[ "; +// for (auto i = 0u; i < M; ++i) +// std::for_each(begin(a), end(a), []() { +// o << std::setw(10) << std::setprecision(3) << std::setfill(' ') << a[i]; +//}); +// o << " ]" << std::endl; + +// return o; +//} + +#endif diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/newtonsMethod/Makefile b/MATH5620_NumericalSolutionsOfDifferentialEquations/newtonsMethod/Makefile new file mode 100644 index 0000000..06119b2 --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/newtonsMethod/Makefile @@ -0,0 +1,13 @@ +OBJS = ../matrix/matrix.hpp ../matrix/matrix_util.hpp ../matrix/vector_util.hpp main.cpp + +CC = g++ +FLAGS = -std=c++17 +DEBUG_FLAGS = -g3 -Og +RELEASE_FLAGS = -g0 -O3 +WARNINGS = -Wall -Wextra -Werror + +release: $(OBJS) + $(CC) $(FLAGS) $(RELEASE_FLAGS) $(WARNINGS) $(OBJS) -o newtonsMethod.out + +debug: $(OBJS) + $(CC) $(FLAGS) $(DEBUG_FLAGS) $(WARNINGS) $(OBJS) -o debug.out diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/newtonsMethod/main.cpp b/MATH5620_NumericalSolutionsOfDifferentialEquations/newtonsMethod/main.cpp new file mode 100644 index 0000000..e97758e --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/newtonsMethod/main.cpp @@ -0,0 +1,6 @@ +#include "explicitEuler.hpp" +#include + +int main() { + std::cout << explicit_euler(0.0, -1.0, 0.4, 0.1, [](double a, double b){return a*a+2*b;}) << std::endl; +} diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/newtonsMethod/manual_newtons_method.md b/MATH5620_NumericalSolutionsOfDifferentialEquations/newtonsMethod/manual_newtons_method.md new file mode 100644 index 0000000..dabc760 --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/newtonsMethod/manual_newtons_method.md @@ -0,0 +1,49 @@ +--- +title: Newton's Method +math: true +layout: default +--- + +{% include mathjax.html %} + + Table of Contents +# Newton's Method + +**Routine Name:** `newtons_method` + +**Author:** Philip Nelson + +**Language:** C++ + +## Description + +`newtons_method`, In numerical analysis, Newton's method (also known as the Newton–Raphson method), named after Isaac Newton and Joseph Raphson, is a method for finding successively better approximations to the roots (or zeroes) of a real-valued function. It is one example of a root-finding algorithm. [1](https://en.wikipedia.org/wiki/Newton%27s_method) + +## Input + +`newtons_method(T x0, T y0, T x, T dt, F f)` requires: + +* `F f` - the function +* `T dt` - the timestep +* `T x0` - the initial guess +* `const uint MAX_ITERATONS` - the number of iterations + +## Output + +The zero of the function at `x`. + +## Code +{% highlight c++ %} +template +T newtons_method(F f, T dt, T x0, const unsigned int MAX_ITERATONS = 100) +{ + auto tol = maceps().maceps, i = 0u; + while (std::abs(f(x0) - 0) > tol && ++i < MAX_ITERATONS) + { + x0 = x0 - f(x0) / ((f(x0 + dt) - f(x0)) / dt); + } + return x0; +} +{% endhighlight %} + +**Last Modification date:** 3 April 2018 diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/newtonsMethod/newtonsMethod.hpp b/MATH5620_NumericalSolutionsOfDifferentialEquations/newtonsMethod/newtonsMethod.hpp new file mode 100644 index 0000000..e7714e4 --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/newtonsMethod/newtonsMethod.hpp @@ -0,0 +1,19 @@ +#ifndef NEWTONS_METHOD +#define NEWTONS_METHOD + +#include "../machineEpsilon/maceps.hpp" +#include +#include + +template +T newtons_method(F f, T dt, T x0, const unsigned int MAX_ITERATONS = 100) +{ + auto tol = maceps().maceps, i = 0u; + while (std::abs(f(x0) - 0) > tol && ++i < MAX_ITERATONS) + { + x0 = x0 - f(x0) / ((f(x0 + dt) - f(x0)) / dt); + } + return x0; +} + +#endif diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/predictorCorrector/Makefile b/MATH5620_NumericalSolutionsOfDifferentialEquations/predictorCorrector/Makefile new file mode 100644 index 0000000..a6badfd --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/predictorCorrector/Makefile @@ -0,0 +1,13 @@ +OBJS = ../matrix/matrix.hpp ../matrix/matrix_util.hpp ../matrix/vector_util.hpp main.cpp + +CC = g++ +FLAGS = -std=c++17 +DEBUG_FLAGS = -g3 -Og +RELEASE_FLAGS = -g0 -O3 +WARNINGS = -Wall -Wextra -Werror + +release: $(OBJS) + $(CC) $(FLAGS) $(RELEASE_FLAGS) $(WARNINGS) $(OBJS) -o predictorCorrector.out + +debug: $(OBJS) + $(CC) $(FLAGS) $(DEBUG_FLAGS) $(WARNINGS) $(OBJS) -o debug.out diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/predictorCorrector/main.cpp b/MATH5620_NumericalSolutionsOfDifferentialEquations/predictorCorrector/main.cpp new file mode 100644 index 0000000..e97758e --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/predictorCorrector/main.cpp @@ -0,0 +1,6 @@ +#include "explicitEuler.hpp" +#include + +int main() { + std::cout << explicit_euler(0.0, -1.0, 0.4, 0.1, [](double a, double b){return a*a+2*b;}) << std::endl; +} diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/predictorCorrector/manual_predictor_corrector.md b/MATH5620_NumericalSolutionsOfDifferentialEquations/predictorCorrector/manual_predictor_corrector.md new file mode 100644 index 0000000..e1da5fe --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/predictorCorrector/manual_predictor_corrector.md @@ -0,0 +1,160 @@ +--- +title: Predictor Corrector +math: true +layout: default +--- + +{% include mathjax.html %} + + Table of Contents +# Predictor Corrector + +**Routine Name:** `predictor_corrector` + +**Author:** Philip Nelson + +**Language:** C++ + +## Description + +`predictor_corrector` belong to a class of algorithms designed to integrate ordinary differential equations – to find an unknown function that satisfies a given differential equation. All such algorithms proceed in two steps: + +1. The initial, "prediction" step, starts from a function fitted to the function-values and derivative-values at a preceding set of points to extrapolate ("anticipate") this function's value at a subsequent, new point. + +2. The next, "corrector" step refines the initial approximation by using the predicted value of the function and another method to interpolate that unknown function's value at the same subsequent point. +[1](https://en.wikipedia.org/wiki/Predictor–corrector_method) + +Three families of linear multistep methods are commonly used: Adams–Bashforth methods, Adams–Moulton methods, and the backward differentiation formulas (BDFs). + +The Adams–Bashforth methods are explicit methods. The coefficients are \\(a_{s-1}=-1\\) and \\(a_{s-2}=\cdots =a_0=0\\), while the \\(b_j\\) are chosen such that the methods has order s (this determines the methods uniquely). + +The Adams–Moulton methods are similar to the Adams–Bashforth methods in that they also have \\(a_{s-1}=-1\\) and \\(a_{s-2}=\cdots =a_0=0\\). Again the b coefficients are chosen to obtain the highest order possible. However, the Adams–Moulton methods are implicit methods. By removing the restriction that \\(b_s = 0\\) , an s-step Adams–Moulton method can reach order s+1, while an s-step Adams–Bashforth methods has only order s. +[2](https://en.wikipedia.org/wiki/Linear_multistep_method) + +## Input + +`predictor_corrector(T x0, T y0, T x, T dt, F f)` requires: + +* `T x0` - the initial `x` +* `T y0` - the initial `y` +* `T x` - the value of `x` for which you want to find value of `y` +* `T dt` - the delta t step +* `F f` - the function defining \\(\frac{dy}{dx}\\) + +## Output + +The value of `y` at `x`. + +## Code +{% highlight c++ %} +template +T predictor_corrector(T x0, T y0, T x, T dt, F f, unsigned int n, const unsigned int MAX_ITERATIONS = 1000) +{ + double A, B, ALPHA, H, t, W[MAX_ITERATIONS], K1, K2, K3, K4; + + A = 0.0; + B = 2.0; + + std::cout.setf(std::ios::fixed, std::ios::floatfield); + std::cout.precision(10); + + H = (B - A) / n; + t = A; + W[0] = 0.5; // initial value + + for (auto i = 1u; i <= 3u; i++) + { + K1 = H * f(t, W[i - 1]); + K2 = H * f(t + H / 2.0, W[i - 1] + K1 / 2.0); + K3 = H * f(t + H / 2.0, W[i - 1] + K2 / 2.0); + K4 = H * f(t + H, W[i - 1] + K3); + + W[i] = W[i - 1] + 1 / 6.0 * (K1 + 2.0 * K2 + 2.0 * K3 + K4); + + t = A + i * H; + + std::cout << "At time " << t << " the solution = " << W[i] << std::endl; + } + + for (auto i = 4u; i <= n; i++) + { + K1 = 55.0 * f(t, W[i - 1]) - 59.0 * f(t - H, W[i - 2]) + + 37.0 * f(t - 2.0 * H, W[i - 3]) - 9.0 * f(t - 3.0 * H, W[i - 4]); + W[i] = W[i - 1] + H / 24.0 * K1; + + t = A + i * H; + + std::cout << "At time " << t << " the solution = " << W[i] << std::endl; + } + + return W[i]; +} +{% endhighlight %} + +## Example +{% highlight c++ %} +int main() +{ + std::cout << + predictor_corrector(0.0, -1.0, 0.4, 0.1, [](double a, double b){return a*a+2*b;}) + << std::endl; +} +{% endhighlight %} + +## Result +``` +----- Lambda Differential Equation ----- + +lambda = 1 +exact 0 -> 10 +approx 0 -> 10 +exact 0.2 -> 12.214 +approx 0.2 -> 12.214 +exact 0.4 -> 14.9182 +approx 0.4 -> 14.9183 +exact 0.6 -> 18.2212 +approx 0.6 -> 18.2212 + +lambda = -1 +exact 0 -> 10 +approx 0 -> 10 +exact 0.2 -> 8.18731 +approx 0.2 -> 8.18732 +exact 0.4 -> 6.7032 +approx 0.4 -> 6.70321 +exact 0.6 -> 5.48812 +approx 0.6 -> 5.48813 + +lambda = 100 +exact 0 -> 10 +approx 0 -> 10 +exact 0.2 -> 4.85165e+09 +approx 0.2 -> 4.90044e+09 +exact 0.4 -> 2.35385e+18 +approx 0.4 -> inf +exact 0.6 -> 1.14201e+27 + +----- Logistic Differential Equaiton ----- + +p0 = 25 +exact 0 -> 25 +approx 0 -> 25 +exact 0.2 -> 25.4922 +approx 0.2 -> 25.4922 +exact 0.4 -> 25.9937 +approx 0.4 -> 25.9937 +exact 0.6 -> 26.5049 +approx 0.6 -> 26.5049 + +p0 = 40000 +exact 0 -> 40000 +approx 0 -> 40000 +exact 0.2 -> 22570.2 +approx 0.2 -> 22570.4 +exact 0.4 -> 15815.2 +approx 0.4 -> 15815.4 +exact 0.6 -> 12228 +approx 0.6 -> 12228.2 +``` + +**Last Modification date:** 3 April 2018 diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/predictorCorrector/predictorCorrector.hpp b/MATH5620_NumericalSolutionsOfDifferentialEquations/predictorCorrector/predictorCorrector.hpp new file mode 100644 index 0000000..a9ae001 --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/predictorCorrector/predictorCorrector.hpp @@ -0,0 +1,51 @@ +#ifndef EXPLICIT_EULER_HPP +#define EXPLICIT_EULER_HPP + +#include "../machineEpsilon/maceps.hpp" +#include +#include + +template +T predictor_corrector(T x0, T y0, T x, T dt, F f, unsigned int n, const unsigned int MAX_ITERATIONS = 1000) +{ + double A, B, ALPHA, H, t, W[MAX_ITERATIONS], K1, K2, K3, K4; + + A = 0.0; + B = 2.0; + + std::cout.setf(std::ios::fixed, std::ios::floatfield); + std::cout.precision(10); + + H = (B - A) / n; + t = A; + W[0] = 0.5; // initial value + + for (auto i = 1u; i <= 3u; i++) + { + K1 = H * f(t, W[i - 1]); + K2 = H * f(t + H / 2.0, W[i - 1] + K1 / 2.0); + K3 = H * f(t + H / 2.0, W[i - 1] + K2 / 2.0); + K4 = H * f(t + H, W[i - 1] + K3); + + W[i] = W[i - 1] + 1 / 6.0 * (K1 + 2.0 * K2 + 2.0 * K3 + K4); + + t = A + i * H; + + std::cout << "At time " << t << " the solution = " << W[i] << std::endl; + } + + for (auto i = 4u; i <= n; i++) + { + K1 = 55.0 * f(t, W[i - 1]) - 59.0 * f(t - H, W[i - 2]) + + 37.0 * f(t - 2.0 * H, W[i - 3]) - 9.0 * f(t - 3.0 * H, W[i - 4]); + W[i] = W[i - 1] + H / 24.0 * K1; + + t = A + i * H; + + std::cout << "At time " << t << " the solution = " << W[i] << std::endl; + } + + return W[i]; +} + +#endif diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/rungeKuttaOrder2/Makefile b/MATH5620_NumericalSolutionsOfDifferentialEquations/rungeKuttaOrder2/Makefile new file mode 100644 index 0000000..7d6ab34 --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/rungeKuttaOrder2/Makefile @@ -0,0 +1,13 @@ +OBJS = ../matrix/matrix.hpp ../matrix/matrix_util.hpp ../matrix/vector_util.hpp main.cpp + +CC = g++ +FLAGS = -std=c++17 +DEBUG_FLAGS = -g3 -Og +RELEASE_FLAGS = -g0 -O3 +WARNINGS = -Wall -Wextra -Werror + +release: $(OBJS) + $(CC) $(FLAGS) $(RELEASE_FLAGS) $(WARNINGS) $(OBJS) -o rungeKuttaOrder2.out + +debug: $(OBJS) + $(CC) $(FLAGS) $(DEBUG_FLAGS) $(WARNINGS) $(OBJS) -o debug.out diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/rungeKuttaOrder2/main.cpp b/MATH5620_NumericalSolutionsOfDifferentialEquations/rungeKuttaOrder2/main.cpp new file mode 100644 index 0000000..e97758e --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/rungeKuttaOrder2/main.cpp @@ -0,0 +1,6 @@ +#include "explicitEuler.hpp" +#include + +int main() { + std::cout << explicit_euler(0.0, -1.0, 0.4, 0.1, [](double a, double b){return a*a+2*b;}) << std::endl; +} diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/rungeKuttaOrder2/manual_runge_kutta_order2.md b/MATH5620_NumericalSolutionsOfDifferentialEquations/rungeKuttaOrder2/manual_runge_kutta_order2.md new file mode 100644 index 0000000..8f404c8 --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/rungeKuttaOrder2/manual_runge_kutta_order2.md @@ -0,0 +1,130 @@ +--- +title: Runge Kutta Order 2 +math: true +layout: default +--- + +{% include mathjax.html %} + + Table of Contents +# Kunge Kutta Order 2 + +**Routine Name:** `runge_kutta_order2` + +**Author:** Philip Nelson + +**Language:** C++ + +## Description + +`runge_kutta_order2` is in a family of implicit and explicit iterative methods, which include the well-known routine called the Euler Method, used in temporal discretization for the approximate solutions of ordinary differential equations. These methods were developed around 1900 by the German mathematicians C. Runge and M. W. Kutta.[1](https://en.wikipedia.org/wiki/Runge–Kutta_methods) + +## Input + +`runge_kutta_order2(F f, unsigned int n, T h, T x0, T y0)` requires: + +* `F f` - the function +* `T x0` - the initial `x` +* `T y0` - the initial `y` +* `T h` - the step size + +## Output + +The value of `y` at `x`. + +## Code +{% highlight c++ %} +template +T runge_kutta_order2(F f, unsigned int n, T h, T x0, T y0) +{ + long double y[n], x[n], k[n][n]; + + x[0] = x0; + y[0] = y0; + + for (auto i = 1u; i <= n; i++) + { + x[i] = x[i - 1] + h; + } + + for (auto j = 1u; j <= n; j++) + { + k[1][j] = h * f(x[j - 1], y[j - 1]); + std::cout << "K[1] = " << k[1][j] << "\t"; + k[2][j] = h * f(x[j - 1] + h, y[j - 1] + k[1][j]); + std::cout << "K[2] = " << k[2][j] << "\n"; + y[j] = y[j - 1] + ((k[1][j] + k[2][j]) / 2); + std::cout << "y[" << h * j << "] = " << std::setprecision(5) << y[j] + << std::endl; + } + return y; +} +{% endhighlight %} + +## Example +{% highlight c++ %} +int main() +{ + std::cout << + runge_kutta_order2(0.0, -1.0, 0.4, 0.1, [](double a, double b){return a*a+2*b;}) + << std::endl; +} +{% endhighlight %} + +## Result +``` +----- Lambda Differential Equation ----- + +lambda = 1 +exact 0 -> 10 +approx 0 -> 10 +exact 0.2 -> 12.214 +approx 0.2 -> 12.214 +exact 0.4 -> 14.9182 +approx 0.4 -> 14.9183 +exact 0.6 -> 18.2212 +approx 0.6 -> 18.2212 + +lambda = -1 +exact 0 -> 10 +approx 0 -> 10 +exact 0.2 -> 8.18731 +approx 0.2 -> 8.18732 +exact 0.4 -> 6.7032 +approx 0.4 -> 6.70321 +exact 0.6 -> 5.48812 +approx 0.6 -> 5.48813 + +lambda = 100 +exact 0 -> 10 +approx 0 -> 10 +exact 0.2 -> 4.85165e+09 +approx 0.2 -> 4.90044e+09 +exact 0.4 -> 2.35385e+18 +approx 0.4 -> inf +exact 0.6 -> 1.14201e+27 + +----- Logistic Differential Equaiton ----- + +p0 = 25 +exact 0 -> 25 +approx 0 -> 25 +exact 0.2 -> 25.4922 +approx 0.2 -> 25.4922 +exact 0.4 -> 25.9937 +approx 0.4 -> 25.9937 +exact 0.6 -> 26.5049 +approx 0.6 -> 26.5049 + +p0 = 40000 +exact 0 -> 40000 +approx 0 -> 40000 +exact 0.2 -> 22570.2 +approx 0.2 -> 22570.4 +exact 0.4 -> 15815.2 +approx 0.4 -> 15815.4 +exact 0.6 -> 12228 +approx 0.6 -> 12228.2 +``` + +**Last Modification date:** 3 April 2018 diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/rungeKuttaOrder2/rungeKutta2.hpp b/MATH5620_NumericalSolutionsOfDifferentialEquations/rungeKuttaOrder2/rungeKutta2.hpp new file mode 100644 index 0000000..2d96cde --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/rungeKuttaOrder2/rungeKutta2.hpp @@ -0,0 +1,34 @@ +#ifndef RUNGE_KUTTA_2_HPP +#define RUNGE_KUTTA_2_HPP + +#include "../machineEpsilon/maceps.hpp" +#include +#include + +template +T runge_kutta_order2(F f, unsigned int n, T h, T x0, T y0) +{ + long double y[n], x[n], k[n][n]; + + x[0] = x0; + y[0] = y0; + + for (auto i = 1u; i <= n; i++) + { + x[i] = x[i - 1] + h; + } + + for (auto j = 1u; j <= n; j++) + { + k[1][j] = h * f(x[j - 1], y[j - 1]); + std::cout << "K[1] = " << k[1][j] << "\t"; + k[2][j] = h * f(x[j - 1] + h, y[j - 1] + k[1][j]); + std::cout << "K[2] = " << k[2][j] << "\n"; + y[j] = y[j - 1] + ((k[1][j] + k[2][j]) / 2); + std::cout << "y[" << h * j << "] = " << std::setprecision(5) << y[j] + << std::endl; + } + return x; +} + +#endif diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/rungeKuttaOrder4/Makefile b/MATH5620_NumericalSolutionsOfDifferentialEquations/rungeKuttaOrder4/Makefile new file mode 100644 index 0000000..06119b2 --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/rungeKuttaOrder4/Makefile @@ -0,0 +1,13 @@ +OBJS = ../matrix/matrix.hpp ../matrix/matrix_util.hpp ../matrix/vector_util.hpp main.cpp + +CC = g++ +FLAGS = -std=c++17 +DEBUG_FLAGS = -g3 -Og +RELEASE_FLAGS = -g0 -O3 +WARNINGS = -Wall -Wextra -Werror + +release: $(OBJS) + $(CC) $(FLAGS) $(RELEASE_FLAGS) $(WARNINGS) $(OBJS) -o newtonsMethod.out + +debug: $(OBJS) + $(CC) $(FLAGS) $(DEBUG_FLAGS) $(WARNINGS) $(OBJS) -o debug.out diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/rungeKuttaOrder4/main.cpp b/MATH5620_NumericalSolutionsOfDifferentialEquations/rungeKuttaOrder4/main.cpp new file mode 100644 index 0000000..e97758e --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/rungeKuttaOrder4/main.cpp @@ -0,0 +1,6 @@ +#include "explicitEuler.hpp" +#include + +int main() { + std::cout << explicit_euler(0.0, -1.0, 0.4, 0.1, [](double a, double b){return a*a+2*b;}) << std::endl; +} diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/rungeKuttaOrder4/manual_runge_kutta_order4.md b/MATH5620_NumericalSolutionsOfDifferentialEquations/rungeKuttaOrder4/manual_runge_kutta_order4.md new file mode 100644 index 0000000..3c1ccc4 --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/rungeKuttaOrder4/manual_runge_kutta_order4.md @@ -0,0 +1,130 @@ +--- +title: Runte Kutta Order 4 +math: true +layout: default +--- + +{% include mathjax.html %} + + Table of Contents +# Kunge Kutta Order 4 + +**Routine Name:** `runge_kutta_order4` + +**Author:** Philip Nelson + +**Language:** C++ + +## Description + +`runge_kutta_order4` is in a family of implicit and explicit iterative methods, which include the well-known routine called the Euler Method, used in temporal discretization for the approximate solutions of ordinary differential equations. These methods were developed around 1900 by the German mathematicians C. Runge and M. W. Kutta.[1](https://en.wikipedia.org/wiki/Runge–Kutta_methods) + +## Input + +`runge_kutta_order4(F f, unsigned int n, T h, T x0, T y0)` requires: + +* `F f` - the function +* `T x0` - the initial `x` +* `T y0` - the initial `y` +* `T h` - the step size + +## Output + +The value of `y` at `x`. + +## Code +{% highlight c++ %} +template +T runge_kutta_order4(F f, unsigned int n, T h, T x0, T y0) +{ + long double y[n], x[n], k[n][n]; + + x[0] = x0; + y[0] = y0; + + for (auto i = 1u; i <= n; i++) + { + x[i] = x[i - 1] + h; + } + + for (auto j = 1u; j <= n; j++) + { + k[1][j] = h * f(x[j - 1], y[j - 1]); + std::cout << "K[1] = " << k[1][j] << "\t"; + k[2][j] = h * f(x[j - 1] + h, y[j - 1] + k[1][j]); + std::cout << "K[2] = " << k[2][j] << "\n"; + y[j] = y[j - 1] + ((k[1][j] + k[2][j]) / 2); + std::cout << "y[" << h * j << "] = " << std::setprecision(5) << y[j] + << std::endl; + } + return y; +} +{% endhighlight %} + +## Example +{% highlight c++ %} +int main() +{ + std::cout << + runge_kutta_order4(0.0, -1.0, 0.4, 0.1, [](double a, double b){return a*a+2*b;}) + << std::endl; +} +{% endhighlight %} + +## Result +``` +----- Lambda Differential Equation ----- + +lambda = 1 +exact 0 -> 10 +approx 0 -> 10 +exact 0.2 -> 12.214 +approx 0.2 -> 12.214 +exact 0.4 -> 14.9182 +approx 0.4 -> 14.9183 +exact 0.6 -> 18.2212 +approx 0.6 -> 18.2212 + +lambda = -1 +exact 0 -> 10 +approx 0 -> 10 +exact 0.2 -> 8.18731 +approx 0.2 -> 8.18732 +exact 0.4 -> 6.7032 +approx 0.4 -> 6.70321 +exact 0.6 -> 5.48812 +approx 0.6 -> 5.48813 + +lambda = 100 +exact 0 -> 10 +approx 0 -> 10 +exact 0.2 -> 4.85165e+09 +approx 0.2 -> 4.90044e+09 +exact 0.4 -> 2.35385e+18 +approx 0.4 -> inf +exact 0.6 -> 1.14201e+27 + +----- Logistic Differential Equaiton ----- + +p0 = 25 +exact 0 -> 25 +approx 0 -> 25 +exact 0.2 -> 25.4922 +approx 0.2 -> 25.4922 +exact 0.4 -> 25.9937 +approx 0.4 -> 25.9937 +exact 0.6 -> 26.5049 +approx 0.6 -> 26.5049 + +p0 = 40000 +exact 0 -> 40000 +approx 0 -> 40000 +exact 0.2 -> 22570.2 +approx 0.2 -> 22570.4 +exact 0.4 -> 15815.2 +approx 0.4 -> 15815.4 +exact 0.6 -> 12228 +approx 0.6 -> 12228.2 +``` + +**Last Modification date:** 3 April 2018 diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/rungeKuttaOrder4/rungeKutta4.hpp b/MATH5620_NumericalSolutionsOfDifferentialEquations/rungeKuttaOrder4/rungeKutta4.hpp new file mode 100644 index 0000000..bd336d4 --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/rungeKuttaOrder4/rungeKutta4.hpp @@ -0,0 +1,20 @@ +#ifndef EXPLICIT_EULER_HPP +#define EXPLICIT_EULER_HPP + +#include "../machineEpsilon/maceps.hpp" +#include +#include + +template +T explicit_euler(T x0, T y0, T x, T dt, F f) +{ + auto tol = maceps().maceps; + while (std::abs(x - x0) > tol) + { + y0 = y0 + (dt * f(x0, y0)); + x0 += dt; + } + return y0; +} + +#endif diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/secondOrderLinear/manual.md b/MATH5620_NumericalSolutionsOfDifferentialEquations/secondOrderLinear/manual.md index be99402..0bed91d 100644 --- a/MATH5620_NumericalSolutionsOfDifferentialEquations/secondOrderLinear/manual.md +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/secondOrderLinear/manual.md @@ -6,7 +6,7 @@ layout: default {% include mathjax.html %} - Table of Contents + Table of Contents # Second Order Linear Constant Coefficents **Routine Name:** solcc diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/stabilityAnalysis/laxWendroff.md b/MATH5620_NumericalSolutionsOfDifferentialEquations/stabilityAnalysis/laxWendroff.md new file mode 100644 index 0000000..d2f5bc0 --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/stabilityAnalysis/laxWendroff.md @@ -0,0 +1,22 @@ +--- +title: laxWendroff Analysis +math: true +layout: default +--- + +{% include mathjax.html %} + + Table of Contents +# Lax Wendroff Stability Analysis + +**Author:** Philip Nelson +(see also p.213 of text) + +\\(g(\xi) = 1 - \frac{1}{2}v(e^{i\xi h} - e^{-i\xi h}) + \frac{1}{2}v^2(e^{i\xi h} - 2 + e^{-i\xi h})\\) +\\(= 1 - iv\sin(\xi h) + v^2(\cos(\xi h) - 1)\\) +\\(= 1 - iv \Big( 2\sin(\frac{\xi h}{2})cos(\frac{\xi h}{2} \Big) + v^2 \Big( 2sin^2(\frac{\xi h}{2}) \Big)\\) +\\( \implies |g(\xi)|^2 = \Big( 1 - 2sin^2(\frac{\xi h}{2}) \Big)^2 + 4\sin^2(\frac{\xi h}{2})cos^2(\frac{\xi h}{2})\\) +\\(= 1 - 4v^2(1 - v^2)sin^4(\frac{\xi h}{2})\\) +\\( \implies \text{stable for } |v| \leq 1\\) + +**Last Modification date:** 5 May 2018 diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/stabilityAnalysis/warmingAndBeam.md b/MATH5620_NumericalSolutionsOfDifferentialEquations/stabilityAnalysis/warmingAndBeam.md new file mode 100644 index 0000000..32c89f8 --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/stabilityAnalysis/warmingAndBeam.md @@ -0,0 +1,25 @@ +--- +title: Warming and Beam Analysis +math: true +layout: default +--- + +{% include mathjax.html %} + + Table of Contents +# Warming and Beam Stability Analysis + +**Author:** Philip Nelson +(see also p.213 of text) + +\\(g(\xi) = 1 - \frac{v}{2}(3 - 4e^{-i\xi h} + e^{-i\xi 2h}) + \frac{v^2}{2}(1 - 2e^{-i\xi h} + e^{-i\xi 2h})\\) +\\(\implies e^{i\xi h} g(\xi) = e^{i\xi h} - \frac{v}{2}(3e^{i\xi h} - 4 + e^{-i\xi h}) + \frac{v^2}{2}(e^{i\xi h} - 2 + e^{-i\xi h})\\) +\\(= e^{i\xi h} - \frac{v}{2}(3e^{i\xi h} - 4 + e^{-i\xi h}) - v^2 + v^2\cos(\xi h)\\) +\\(= e^{i\xi h} - \frac{v}{2}(2e^{i\xi h} - 4) + v^2(\cos(\xi h) - 1)\\) +\\(= e^{i\xi h}(1 - v) - 2v + v\cos(\xi h) + v^2(\cos(\xi h) - 1)\\) +\\(= e^{i\xi h}(1 - v) - v(2 - \cos(\xi h)) + v^2(\cos(\xi h) - 1)\\) +\\(= (\cos(\xi h) + i\sin(\xi h))(1 - v) - v(2 - \cos(\xi h)) + v^2(\cos(\xi h) - 1)\\) +\\(\implies |g(\xi)|^2 = sin^2(\xi h) + \Big( \cos(\xi h)(1 - v) - v(2 - \cos(\xi h)) + v^2(\cos(\xi h) - 1) \Big)^2\\) +\\(\implies \text{stable for } 0 \leq v \leq 2\\) + +**Last Modification date:** 5 May 2018 diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/template.md b/MATH5620_NumericalSolutionsOfDifferentialEquations/template.md new file mode 100644 index 0000000..40959c3 --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/template.md @@ -0,0 +1,54 @@ +--- +title: <++> +math: true +layout: default +--- + +{% include mathjax.html %} + + Table of Contents +# <++> title + +**Routine Name:** <++> + +**Author:** Philip Nelson + +**Language:** C++ + +## Description + +<++> +`pNorm` calculates the p norm of a vector + +## Input + +<++> +`pNorm(std::array v, unsigned int const& p)` requires: + +* `std::array v` - a column vector of type `T` and size `M` +* `unsigned int p` - the p norm that is desired + +## Output + +<++> +A double with the desired norm value. + +## Code +{% highlight c++ %} +<++> +{% endhighlight %} + +## Example +{% highlight c++ %} +int main() +{ +<++> +} +{% endhighlight %} + +## Result +``` +<++> +``` + +**Last Modification date:** <++>10 February 2018 diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/testConjugateGradientFivePoint/Makefile b/MATH5620_NumericalSolutionsOfDifferentialEquations/testConjugateGradientFivePoint/Makefile new file mode 100644 index 0000000..b5da4ee --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/testConjugateGradientFivePoint/Makefile @@ -0,0 +1,12 @@ +OBJS = main.cpp ../matrix/matrix.hpp + +CC = g++ +FLAGS = -std=c++17 +DEBUG_FLAGS = -Og -g3 -Wall -Wextra +RELEASE_FLAGS = -O3 + +release: $(OBJS) + $(CC) $(RELEASE_FLAGS) $(FLAGS) $(OBJS) -o test.out + +debug: $(OBJS) + $(CC) $(DEBUG_FLAGS) $(FLAGS) $(OBJS) -o debug.out diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/testConjugateGradientFivePoint/main.cpp b/MATH5620_NumericalSolutionsOfDifferentialEquations/testConjugateGradientFivePoint/main.cpp new file mode 100644 index 0000000..03b315f --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/testConjugateGradientFivePoint/main.cpp @@ -0,0 +1,11 @@ +#include "../matrix/matrix.hpp" +#include "../matrix/matrix_util.hpp" +#include "../conjugateGradient/conjugateGradient.hpp" +#include + +int main() +{ + auto answer = solveFivePointStencil(0.0, 1.0, sin); + auto finalMat = arrayToMat(answer); + std::cout << "Answer in Matrix Form\n" << finalMat << std::endl; +} diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/testConjugateGradientFivePoint/manual_solve_five_point_stencil_test.md b/MATH5620_NumericalSolutionsOfDifferentialEquations/testConjugateGradientFivePoint/manual_solve_five_point_stencil_test.md new file mode 100644 index 0000000..b8b106e --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/testConjugateGradientFivePoint/manual_solve_five_point_stencil_test.md @@ -0,0 +1,98 @@ +--- +title: 5-Point Stencil Test +math: true +layout: default +--- + +{% include mathjax.html %} + + Table of Contents +# Five Point Stencil Test + +**Routine Name:** solveFivePointStencil + +**Author:** Philip Nelson + +**Language:** C++ + +## Description + +`fivePointStencil` solves Poisson's Equation for an arbitrary mesh size and forcing function + +## Input + +`solveFivePointStencil(T a, T b, F f)` requires: + +* `T a` - the start boundary +* `T b` - the end boundary +* `F f` - the forcing function +* `T` - template param for the type +* `N` - template param for the dimension of the mesh + +## Output + +An array with the approximated solution + +## Code +{% highlight c++ %} +template +auto solveFivePointStencil(T a, T b, F f) +{ + auto mesh = generateMesh(a, b); + auto bv = initMeshB(mesh, f); + auto stencil = fivePointStencil(); + auto res = conjugate_gradient(stencil, bv); + return res; +} +{% endhighlight %} + +solveFivePointStencil relies on: +[fivePointStencil](./manual_gen_five_point_stencil)| +[generateMesh](./manual_gen_mesh)| +[initMeshB](./manual_init_b)| +[conjugate_gradient](../conjugateGradient/manual_conjugate_gradient)| + +## Example +{% highlight c++ %} +int main() +{ + auto answer = solveFivePointStencil(0.0, 1.0, sin); + auto finalMat = arrayToMat(answer); + std::cout << "Answer in Matrix Form\n" << finalMat << std::endl; +} +{% endhighlight %} + +## Result +``` +mesh +| 0.0625 0.125 0.188 | +| 0.125 0.25 0.375 | +| 0.188 0.375 0.562 | + + +stencil +| -4 1 0 1 0 0 0 0 0 | +| 1 -4 1 0 1 0 0 0 0 | +| 0 1 -4 0 0 1 0 0 0 | +| 1 0 0 -4 1 0 1 0 0 | +| 0 1 0 1 -4 1 0 1 0 | +| 0 0 1 0 1 -4 0 0 1 | +| 0 0 0 1 0 0 -4 1 0 | +| 0 0 0 0 1 0 1 -4 1 | +| 0 0 0 0 0 1 0 1 -4 | + + +bv +[ 0.0625 0.125 0.186 0.125 0.247 0.366 0.186 0.366 0.533 ] + + +result +[ -0.097 -0.163 -0.154 -0.163 -0.276 -0.266 -0.154 -0.266 -0.266 ] + +Answer in Matrix Form +| -0.097 -0.163 -0.154 | +| -0.163 -0.276 -0.266 | +| -0.154 -0.266 -0.266 | +``` + +**Last Modification date:** 21 March 2018 diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/upwinding/Makefile b/MATH5620_NumericalSolutionsOfDifferentialEquations/upwinding/Makefile new file mode 100644 index 0000000..f26aed4 --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/upwinding/Makefile @@ -0,0 +1,13 @@ +OBJS = ../matrix/matrix.hpp ../matrix/matrix_util.hpp ../matrix/vector_util.hpp main.cpp + +CC = g++ +FLAGS = -std=c++17 +DEBUG_FLAGS = -g3 -Og +RELEASE_FLAGS = -g0 -O3 +WARNINGS = -Wall -Wextra -Werror + +release: $(OBJS) + $(CC) $(FLAGS) $(RELEASE_FLAGS) $(WARNINGS) $(OBJS) -o upwinding.out + +debug: $(OBJS) + $(CC) $(FLAGS) $(DEBUG_FLAGS) $(WARNINGS) $(OBJS) -o debug.out diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/upwinding/main.cpp b/MATH5620_NumericalSolutionsOfDifferentialEquations/upwinding/main.cpp new file mode 100644 index 0000000..43298bf --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/upwinding/main.cpp @@ -0,0 +1,32 @@ +#include "upwinding.hpp" +#include +#include + +int main() +{ + constexpr double xDomain[] = {0.0, 1.0}; + constexpr double tDomain[] = {0.0, 1.0}; + // constexpr auto dt = 1.0e-3; + // constexpr auto dx = 1.0e-3; + constexpr auto dt = .1; + constexpr auto dx = .1; + constexpr std::size_t S = ((xDomain[1] - xDomain[0]) / dx); + auto c = 0.7; + auto eta = [](const double& x) -> double { + return (x >= 0.3 && x <= 0.6) ? 100 : 0; + }; + + auto solution = upwinding(xDomain, tDomain, dx, dt, eta, c); + + for (auto i = 0u; i < solution.size(); ++i) + { + for (auto j = 0u; j < solution[i].size(); ++j) + { + std::cout << std::setprecision(3) << std::setw(7) << solution[j][i] + << " "; + } + std::cout << '\n'; + } + + return EXIT_SUCCESS; +} diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/upwinding/manual_upwinding.md b/MATH5620_NumericalSolutionsOfDifferentialEquations/upwinding/manual_upwinding.md new file mode 100644 index 0000000..e319103 --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/upwinding/manual_upwinding.md @@ -0,0 +1,127 @@ +--- +title: Upwinding +math: true +layout: default +--- + +{% include mathjax.html %} + + Table of Contents +# Upwinding + +**Routine Name:** `upwinding` + +**Author:** Philip Nelson + +**Language:** C++ + +## Description + +`upwinding` is in a class of numerical discretization methods for solving hyperbolic partial differential equations. Upwind schemes use an adaptive or solution-sensitive finite difference stencil to numerically simulate the direction of propagation of information in a flow field. The upwind schemes attempt to discretize hyperbolic partial differential equations by using differencing biased in the direction determined by the sign of the characteristic speeds. Historically, the origin of upwind methods can be traced back to the work of Courant, Isaacson, and Rees who proposed the CIR method.[1](https://en.wikipedia.org/wiki/Upwind_scheme) + +## Input + +{% highlight c++ %} +upwinding(const T xDomain[], + const T tDomain[], + const T dx, + const T dt, + F eta, + const T c + ) +{% endhighlight %} + +* `T xDomain[]` - two member array with the spacial bounds +* `T tDomain[]` - two member array with the temporal bounds +* `T dx` - the spacial step +* `T dt` - the temporal step +* `F eta` - the function eta +* `T c` - the constant + +## Output + +A `std::vector>` with the solution over time. + +## Code +{% highlight c++ %} +template +std::vector> upwinding(const T xDomain[], + const T tDomain[], + const T dx, + const T dt, + F eta, + const T c) +{ + std::array U; + for (auto i = 0u; i < S - 2; ++i) + { + auto A = xDomain[0]; + auto x = A + (i + 1) * dx; + U[i] = eta(x); + } + + std::vector coeffs = {c * dx / dt, 1 - c * dx / dt}; + + auto A = makeNDiag(coeffs, -1u); + + std::vector> solution = {U}; + + for (auto t = tDomain[0]; t < tDomain[1]; t += dt) + { + U = A * U; + solution.push_back(U); + } + + return solution; +} +{% endhighlight %} + +## Example +{% highlight c++ %} +int main() +{ + constexpr double xDomain[] = {0.0, 1.0}; + constexpr double tDomain[] = {0.0, 1.0}; + // constexpr auto dt = 1.0e-3; + // constexpr auto dx = 1.0e-3; + constexpr auto dt = .1; // for view able output + constexpr auto dx = .1; // for view able output + constexpr std::size_t S = ((xDomain[1] - xDomain[0]) / dx); + auto c = 0.7; + auto eta = [](const double& x) -> double { + return (x >= 0.3 && x <= 0.6) ? 100 : 0; + }; + + auto solution = upwinding(xDomain, tDomain, dx, dt, eta, c); + + for (auto i = 0u; i < solution.size(); ++i) + { + for (auto j = 0u; j < solution[i].size(); ++j) + { + std::cout << std::setprecision(3) << std::setw(7) << solution[j][i] + } + std::cout << '\n'; + } + + return EXIT_SUCCESS; +} +{% endhighlight %} + +## Result +``` + 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 + 100 30 9 2.7 0.81 0.243 0.0729 0.0219 + 100 100 51 21.6 8.37 3.08 1.09 0.379 + 100 100 100 65.7 34.8 16.3 7.05 2.88 + 0 70 91 97.3 75.2 46.9 25.5 12.6 + 0 0 49 78.4 91.6 80.1 56.9 34.9 + 0 0 0 34.3 65.2 83.7 81.2 64.2 + 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 + 30 9 2.7 0.81 0.243 0.0729 0.0219 0.00656 + 100 51 21.6 8.37 3.08 1.09 0.379 0.129 + +``` + +**Last Modification date:** 5 May 2018 diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/upwinding/upwinding.hpp b/MATH5620_NumericalSolutionsOfDifferentialEquations/upwinding/upwinding.hpp new file mode 100644 index 0000000..723d780 --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/upwinding/upwinding.hpp @@ -0,0 +1,40 @@ +#ifndef UPWINDING_HPP +#define UPWINDING_HPP + +#include "../matrix/matrix.hpp" +#include "../matrix/matrix_util.hpp" +#include +#include + +template +std::vector> upwinding(const T xDomain[], + const T tDomain[], + const T dx, + const T dt, + F eta, + const T c) +{ + std::array U; + for (auto i = 0u; i < S - 2; ++i) + { + auto A = xDomain[0]; + auto x = A + (i + 1) * dx; + U[i] = eta(x); + } + + std::vector coeffs = {c * dx / dt, 1 - c * dx / dt}; + + auto A = makeNDiag(coeffs, -1u); + + std::vector> solution = {U}; + + for (auto t = tDomain[0]; t < tDomain[1]; t += dt) + { + U = A * U; + solution.push_back(U); + } + + return solution; +} + +#endif diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/warmingAndBeam/Makefile b/MATH5620_NumericalSolutionsOfDifferentialEquations/warmingAndBeam/Makefile new file mode 100644 index 0000000..d7ec507 --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/warmingAndBeam/Makefile @@ -0,0 +1,13 @@ +OBJS = ../matrix/matrix.hpp ../matrix/matrix_util.hpp ../matrix/vector_util.hpp main.cpp + +CC = g++ +FLAGS = -std=c++17 +DEBUG_FLAGS = -g3 -Og +RELEASE_FLAGS = -g0 -O3 +WARNINGS = -Wall -Wextra -Werror + +release: $(OBJS) + $(CC) $(FLAGS) $(RELEASE_FLAGS) $(WARNINGS) $(OBJS) -o warmingAndBeam.out + +debug: $(OBJS) + $(CC) $(FLAGS) $(DEBUG_FLAGS) $(WARNINGS) $(OBJS) -o debug.out diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/warmingAndBeam/main.cpp b/MATH5620_NumericalSolutionsOfDifferentialEquations/warmingAndBeam/main.cpp new file mode 100644 index 0000000..f66ced8 --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/warmingAndBeam/main.cpp @@ -0,0 +1,32 @@ +#include "warmingAndBeam.hpp" +#include +#include + +int main() +{ + constexpr double xDomain[] = {0.0, 1.0}; + constexpr double tDomain[] = {0.0, 1.0}; + // constexpr auto dt = 1.0e-3; + // constexpr auto dx = 1.0e-3; + constexpr auto dt = .1; + constexpr auto dx = .1; + constexpr std::size_t S = ((xDomain[1] - xDomain[0]) / dx); + auto c = 0.7; + auto eta = [](const double& x) -> double { + return (x >= 0.3 && x <= 0.6) ? 100 : 0; + }; + + auto solution = warming_and_beam(xDomain, tDomain, dx, dt, eta, c); + + for (auto i = 0u; i < solution.size(); ++i) + { + for (auto j = 0u; j < solution[i].size(); ++j) + { + std::cout << std::setprecision(3) << std::setw(10) << solution[j][i] + << " "; + } + std::cout << '\n'; + } + + return EXIT_SUCCESS; +} diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/warmingAndBeam/manual_warming_and_beam.md b/MATH5620_NumericalSolutionsOfDifferentialEquations/warmingAndBeam/manual_warming_and_beam.md new file mode 100644 index 0000000..fbae561 --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/warmingAndBeam/manual_warming_and_beam.md @@ -0,0 +1,129 @@ +--- +title: Warming and Beam +math: true +layout: default +--- + +{% include mathjax.html %} + + Table of Contents +# Warming and Beam + +**Routine Name:** `warming_and_beam` + +**Author:** Philip Nelson + +**Language:** C++ + +## Description + +`warming_and_beam` introduced in 1978 by Richard M. Beam and R. F. Warming, is a second order accurate implicit scheme, mainly used for solving non-linear hyperbolic equation. **It is not used much nowadays**.[1](https://en.wikipedia.org/wiki/Lax–Wendroff_method) + +## Input + +{% highlight c++ %} +warming_and_beam(const T xDomain[], + const T tDomain[], + const T dx, + const T dt, + F eta, + const T c + ) +{% endhighlight %} + +* `T xDomain[]` - two member array with the spacial bounds +* `T tDomain[]` - two member array with the temporal bounds +* `T dx` - the spacial step +* `T dt` - the temporal step +* `F eta` - the function eta +* `T c` - the constant + +## Output + +A `std::vector>` with the solution over time. + +## Code +{% highlight c++ %} +template +std::vector> warming_and_beam(const T xDomain[], + const T tDomain[], + const T dx, + const T dt, + F eta, + const T c) +{ + std::array U; + for (auto i = 0u; i < S - 2; ++i) + { + auto A = xDomain[0]; + auto x = A + (i + 1) * dx; + U[i] = eta(x); + } + + auto a = c * dx / dt / dt; + auto b = c * c * dx * dx / dt / dt / 2; + std::vector coeffs = {b - a, 4 * a - 2 * b, 1 - 2 * a + b}; + + auto A = makeNDiag(coeffs, -2u); + + std::vector> solution = {U}; + + for (auto t = tDomain[0]; t < tDomain[1]; t += dt) + { + U = A * U; + solution.push_back(U); + } + + return solution; +} +{% endhighlight %} + +## Example +{% highlight c++ %} +int main() +{ + constexpr double xDomain[] = {0.0, 1.0}; + constexpr double tDomain[] = {0.0, 1.0}; + // constexpr auto dt = 1.0e-3; + // constexpr auto dx = 1.0e-3; + constexpr auto dt = .1; + constexpr auto dx = .1; + constexpr std::size_t S = ((xDomain[1] - xDomain[0]) / dx); + auto c = 0.7; + auto eta = [](const double& x) -> double { + return (x >= 0.3 && x <= 0.6) ? 100 : 0; + }; + + auto solution = warming_and_beam(xDomain, tDomain, dx, dt, eta, c); + + for (auto i = 0u; i < solution.size(); ++i) + { + for (auto j = 0u; j < solution[i].size(); ++j) + { + std::cout << std::setprecision(3) << std::setw(10) << solution[j][i] + << " "; + } + std::cout << '\n'; + } + + return EXIT_SUCCESS; +} +{% endhighlight %} + +## Result +``` + 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 + 100 -1.28e+03 1.63e+04 -2.08e+05 2.65e+06 -3.38e+07 4.31e+08 -5.49e+09 + 100 1.48e+03 -5.39e+04 1.14e+06 -2.02e+07 3.3e+08 -5.14e+09 7.74e+10 + 100 800 3.9e+04 -2.09e+06 5.93e+07 -1.33e+09 2.63e+10 -4.79e+11 + 0 2.08e+03 -1.44e+04 1.62e+06 -8.59e+07 2.86e+09 -7.53e+10 1.72e+12 + 0 -675 6.03e+04 -1.43e+06 7.7e+07 -3.74e+09 1.35e+11 -3.98e+12 + 0 0 -3.26e+04 2.17e+06 -7.8e+07 3.69e+09 -1.69e+11 6.4e+12 + 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 + -1.28e+03 1.63e+04 -2.08e+05 2.65e+06 -3.38e+07 4.31e+08 -5.49e+09 7.01e+10 + 1.48e+03 -5.39e+04 1.14e+06 -2.02e+07 3.3e+08 -5.14e+09 7.74e+10 -1.14e+12 +``` + +**Last Modification date:** 5 May 2018 diff --git a/MATH5620_NumericalSolutionsOfDifferentialEquations/warmingAndBeam/warmingAndBeam.hpp b/MATH5620_NumericalSolutionsOfDifferentialEquations/warmingAndBeam/warmingAndBeam.hpp new file mode 100644 index 0000000..461b941 --- /dev/null +++ b/MATH5620_NumericalSolutionsOfDifferentialEquations/warmingAndBeam/warmingAndBeam.hpp @@ -0,0 +1,42 @@ +#ifndef LAX_WENDORF_HPP +#define LAX_WENDORF_HPP + +#include "../matrix/matrix.hpp" +#include "../matrix/matrix_util.hpp" +#include +#include + +template +std::vector> warming_and_beam(const T xDomain[], + const T tDomain[], + const T dx, + const T dt, + F eta, + const T c) +{ + std::array U; + for (auto i = 0u; i < S - 2; ++i) + { + auto A = xDomain[0]; + auto x = A + (i + 1) * dx; + U[i] = eta(x); + } + + auto a = c * dx / dt / dt; + auto b = c * c * dx * dx / dt / dt / 2; + std::vector coeffs = {b - a, 4 * a - 2 * b, 1 - 2 * a + b}; + + auto A = makeNDiag(coeffs, -2u); + + std::vector> solution = {U}; + + for (auto t = tDomain[0]; t < tDomain[1]; t += dt) + { + U = A * U; + solution.push_back(U); + } + + return solution; +} + +#endif