computational-physics/project1.cpp at master · nhofield/computational-physics · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
#include <iostream>
#include <cstdlib>
#include <cmath>
#include <fstream>
#include <chrono>
#include <ratio>
#include <ctime>
#include <algorithm>

using namespace std;
inline double f(double x){return 100.0*exp(-10.0*x);}
inline double analytical(double x){return 1.0-(1-exp(-10.0))*x - exp(-10.0*x);}

int main(int argc, char* argv[])
  {

    int n = atoi(argv[1]); //input number of grid points

    //Allocating memory for 6 vectors; a, b, c, b_tilde, d_tilde and v (solution)
    double* a       = new double[n-1];
    double* d       = new double[n];
    double* c       = new double[n-1];

    double* u       = new double[n+2];  //Analytical solution
    double* b_tilde = new double[n];
    double* d_tilde = new double[n];
    double* v       = new double[n+2];
    double* b       = new double[n];

    double h = (double)1/(n+1);
    double hh = h*h;

    //Filling a and c
    for(int i = 0; i < n-1; i++)
      {
        a[i] = c[i] = -1.0;
      }
    //Filling d and b
    for(int i = 0; i < n; i++)
      {
        d[i] = 2.0;
        b[i] = hh*f(h*(i+1));
      }
    //Computing analytical solution
    for(int i = 1; i <= n; i++)
      {
        u[i] = analytical(i*h);
      }

    //Boundary conditions
    u[0] = u[n+1] = 0.0;
    v[0] = v[n+1] = 0.0;

    //b_tilde and d_tilde initial conditions
    b_tilde[0] = b[0];
    d_tilde[0] = d[0];

    /*START Algorithm*/

    auto start = chrono::high_resolution_clock::now();

    /*Forward substitution*/
    for(int i = 1; i < n; i++)
      {
        d_tilde[i] = d[i] - (a[i-1]*c[i-1])/d_tilde[i-1];
        b_tilde[i] = b[i] - (a[i-1]*b_tilde[i-1])/d_tilde[i-1];
      }

    //Compute v_n
    v[n] = b_tilde[n-1]/d_tilde[n-1];

      /*Backward substitution*/
    for(int i = n-1 ; i >= 1; i--)
      {
        v[i] = (b_tilde[i-1] - c[i-1]*v[i+1])/d_tilde[i-1];
      }

    /*STOP Algorithm*/

    auto stop = chrono::high_resolution_clock::now();
    auto diff = stop-start;
    cout  << " Runtime of algorithm = " << chrono::duration <double,milli> (diff).count() << "ms" << endl;

    //Not used anymore
    delete []c;
    delete []b_tilde;
    delete []d_tilde;
    delete []a;
    delete []d;
    delete []b;
    //Calc eps_relative

    //Writing result to file
    ofstream outfile;
    outfile.open("proj1data.dat");

    if(*argv[2] == 'Y')
      {
        for(int i = 1; i <= n; i++) outfile << i*h << " " << u[i] << " " << v[i] << endl;
      }
    else
      {
        //Alocating memory for storing relative error
        double* eps_relative = new double[n];
        for(int i = 0; i <= n; i++) eps_relative[i] = log10( abs( ( v[i+1] - u[i+1] )/u[i+1] ) );

        //Fetching maximum value in array, then deallocating memory
        outfile << *max_element(eps_relative, eps_relative + n) << endl;
        delete []eps_relative;
      }
    outfile.close();

    return 0;
  }