Double Sweep and Thomas Algorithm generate different result

M

Michael CQF

Joined
5/28/23
Messages
14
Points
13
I'm learning from book

Financial Instrument Pricing Using C++, 2nd Edition​


Below are my sample test code with real examble from option pricing

probem:
Thomans Algorithm yield different result against double sweep when BCR is not 0.

I have a guess it comes from boundary handling, also have a question, where is a[0] and c[Last one] getting use in Thomas Algorithm?

C++: 
#include <iostream>
#include <vector>
#include "LUSolver.hpp"
#include "DoubleSweep.hpp"
#include "print_utils.hpp"

int main(){
    //std::size_t J = 4;
 
    //double h = 1.0 / static_cast<double>(J);
    // Boundary conditions
    double BCL = 0.0;   // LHS
    //double BCR = 20.19801326693245; // Option real payoff,  but even use a rounding value 20 below will lead to different value
    double BCR = 20.0;
    // Double Sweep
    // std::vector<double> a(J+1,1.0);
    // std::vector<double> b(J+1,-2.0);
    // std::vector<double> c(J+1,1.0);
    // std::vector<double> r(J+1, -2.0*h*h);    // Right-hand side
 
    std::vector<double> a = {0.0, 0.0352, 0.7031, 2.0039, 0.0};
    std::vector<double> b = {0.0, -57.4453, -59.3438, -62.5078, 0.0};
    std::vector<double> c = {0.0, 0.5977, 1.8281, 3.6914, 0.0};
    std::vector<double> r = {0.0, -2.9883, -288.6328, -708.7500, 0.0};
    std::vector<double> a2(a.begin() + 1, a.end() - 1);
    std::vector<double> b2(b.begin() + 1, b.end() - 1);
    std::vector<double> c2(c.begin() + 1, c.end() - 1);
    std::vector<double> r2(r.begin() + 1, r.end() - 1);
    r2[0] -= BCL;
    r2[r2.size() - 1] -= BCR;

 
    LUTridiagonalSolver<double> mySolver3(a2, b2, c2, r2);
    DoubleSweep<double> mySolver(a, b, c, r, BCL, BCR);
    print("solution using Thomas LU Decomposition :{}\n", mySolver3.solve());
    //solution using Thomas LU Decomposition :[0.1064, 5.2293, 11.8262, ]
    print("solution using Double Sweep : {}\n", mySolver.solve());
    //solution using Double Sweep : [0.0, 0.1067, 5.2559, 12.6882, 20.0, ]
 
    return 0;
}
;
 

Attachments

Last edited:
I would first uncouple the 2 algos .. no mixing of a and a2 etc.

Thomas BC taken care of in RHS r term.

J = 4

DS a(J+1)
TA a2(J-1), better than cryptic a2(,,,,).

I use DS in many apps. It originated in former Soviet Union.
It is about 20% faster than TA.
 
Code: 
    std::vector<double> r2(r.begin() + 1, r.end() - 1); => FIRST INDEX is ONE??
    r2[0] -= BCL;
    r2[r2.size() - 1] -= BCR;


edge condition?

r2(J-1)?

take example p. 403 Duffy 2018.
 
Last edited:
Daniel Duffy said:
Code: 
    std::vector<double> r2(r.begin() + 1, r.end() - 1); => FIRST INDEX is ONE??
    r2[0] -= BCL;
    r2[r2.size() - 1] -= BCR;


edge condition?

r2(J-1)?

take example p. 403 Duffy 2018.
Click to expand...

Assuming this is input for Double Sweep

C++: 
   std::vector<double> a = {0.0, 0.0352, 0.7031, 2.0039, 0.0};
    std::vector<double> b = {0.0, -57.4453, -59.3438, -62.5078, 0.0};
    std::vector<double> c = {0.0, 0.5977, 1.8281, 3.6914, 0.0};
    std::vector<double> r = {0.0, -2.9883, -288.6328, -708.7500, 0.0};

What shall be equivalent input of thomas Algorithm? My understanding from textbook is remove the first element of each array and last element of each array

The original textbook use below or double sweep algorithm

Chapter 13 Numerical Linear Algebra\BasicAlgebra\Algebra\TestTridiagonalSolvers.cpp

C++: 
    // Double Sweep
    std::vector<double> a(J+1,1.0);
    std::vector<double> b(J+1,-2.0);
    std::vector<double> c(J+1,1.0);
    std::vector<double> r(J+1, -2.0*h*h);    // Right-hand side


and use below input for Thomas Algorithm

Code: 
    // Thomas algorithm
    std::vector<double> a2(J-1, 1.0);
    std::vector<double> b2(J-1, -2.0);
    std::vector<double> c2(J-1, 1.0);
    std::vector<double> r2(J-1, -2.0*h*h);    // Right-hand side

    // Take the boundary conditions into consideration
    r2[0] -= BCL;
    r2[r2.size()-1] -=  BCR;

However, above is a idea case where all element in a, b, c are same. but what shall I slice if we face a real eample as my one

Code: 
  std::vector<double> a = {0.0, 0.0352, 0.7031, 2.0039, 0.0};

Which is used to solve for crank nicolson
 
I guess the solution is

instead of below from textbook​

C++: 
  r2[0] -= BCL;
   r2[r2.size() - 1] -= BCR;

we adjust like

C++: 
    r2[0] -= BCL*a2[0];

    r2[r2.size() - 1] -= BCR*c2[c2.size()-1];


Full Code is

C++: 
#include <iostream>
#include <vector>
#include "LUSolver.hpp"
#include "DoubleSweep.hpp"
#include "print_utils.hpp"

int main(){
    //std::size_t J = 4;
 
    //double h = 1.0 / static_cast<double>(J);
    // Boundary conditions
    double BCL = 0.0;   // LHS
    //double BCR = 20.19801326693245; // Option real payoff,  but even use a rounding value 20 below will lead to different value
    double BCR = 20.0;
    // Double Sweep
    // std::vector<double> a(J+1,1.0);
    // std::vector<double> b(J+1,-2.0);
    // std::vector<double> c(J+1,1.0);
    // std::vector<double> r(J+1, -2.0*h*h);    // Right-hand side
 
    std::vector<double> a = {0.0, 0.0352, 0.7031, 2.0039, 0.0};
    std::vector<double> b = {0.0, -57.4453, -59.3438, -62.5078, 0.0};
    std::vector<double> c = {0.0, 0.5977, 1.8281, 3.6914, 0.0};
    std::vector<double> r = {0.0, -2.9883, -288.6328, -708.7500, 0.0};
    std::vector<double> a2(a.begin() + 1, a.end() - 1);
    std::vector<double> b2(b.begin() + 1, b.end() - 1);
    std::vector<double> c2(c.begin() + 1, c.end() - 1);
    std::vector<double> r2(r.begin() + 1, r.end() - 1);
    r2[0] -= BCL*a2[0];

    r2[r2.size() - 1] -= BCR*c2[c2.size()-1];



    LUTridiagonalSolver<double> mySolver3(a2, b2, c2, r2);
    DoubleSweep<double> mySolver(a, b, c, r, BCL, BCR);
    print("solution using Thomas LU Decomposition :{}\n", mySolver3.solve());
    //solution using Thomas LU Decomposition :[0.1067, 5.2559, 12.6882, ]
    print("solution using Double Sweep : {}\n", mySolver.solve());
    //solution using Double Sweep : [0.0, 0.1067, 5.2559, 12.6882, 20.0,
 
    return 0;
};

Attached Excel to demo
 

Attachments

Last edited:
Daniel Duffy said:
I can try to debug this code but a better idea is to map the TA in section 13.2.2 (13.11 and 13.12) to code or use the code in my book.

mixing a and a2 is confusing for me.

Tridiagonal matrix algorithm - Wikipedia

en.wikipedia.org en.wikipedia.org
Click to expand...
I'm using the sample in chapter 13


C++: 
   Vector a(J+1,1.0);
    Vector b(J+1,-2.0);
    Vector c(J+1,1.0);
    Vector r(J+1, -2.0*h*h);// Right-hand side
    // Thomas algorithm
    Vector a2(J-1, 1.0);
    Vector b2(J-1, -2.0);
    Vector c2(J-1, 1.0);
    Vector r2(J-1, -2.0*h*h);       // Right-hand side
    // Take the boundary conditions into consideration
    r2[0] -= BCL;
    r2[r2.size()-1] -= BCR;


This is original code

What I mean is, it worked and yield same result because BCL is 0

To accomodate a(double sweep) to a2(thomas)

we need to use below instead

C++: 
  r2[0] -= BCL*a2[0];

  r2[r2.size() - 1] -= BCR*c2[c2.size()-1];
 
Back
Top Bottom