Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- #include <iostream>
- #include <random>
- #include <fstream>
- #include <assert.h>
- #include <math.h>
- #include <iomanip> // For std::fixed and std::setprecision
- int main() {
- double k1b, k1f, k2b, k2f;
- double n0, ne;
- double error;
- double dt;
- double time;
- int steps;
- error = 1.0E-3;
- k1f = 2.24957*1.0E-10; k1b = 4.34413*1.0E-10;
- k2f = 7.10391*1.0E-11; k2b = 1.36891*1.0E-9;
- ne = 1.0E+15; n0 = 1.0E+20;
- dt = 1.0E-8; time = 1.0E-1;
- steps = static_cast<int>(time/dt);
- std::vector<std::array<double, 3>> n(steps);
- //initial condition handling
- n[0][0] = n0; n[0][1] = 0.0; n[0][2] = 0.0;
- for (int i = 1; i < steps; i++) {
- n[i][0] = n[i-1][0] + dt*(-(k1f+k2f)*n[i-1][0] + k1b*n[i-1][1] + k2b*n[i-1][2])*ne;
- n[i][1] = n[i-1][1] + dt*(k1f*n[i-1][0] - k1b*n[i-1][1])*ne;
- n[i][2] = n[i-1][2] + dt*(k2f*n[i-1][0] - k2b*n[i-1][2])*ne;
- std::cout << n[i][0] << " " << n[i][1] << " " << n[i][2] << "\n";
- double rel0, rel1, rel2;
- rel0 = n[i][0]/n[i-1][0]; rel1 = n[i][1]/n[i-1][1]; rel2 = n[i][2]/n[i-1][2];
- if (std::abs(n[i-1][1]) < 1e-30) rel1 = 1.0;
- if (std::abs(n[i-1][2]) < 1e-30) rel2 = 1.0;
- if ((std::abs(rel0) < error) && (std::abs(rel1) < error) && (std::abs(rel2) < error)) {
- std::cout << "Accuracy reached at step: " << i << "\n";
- break;
- }
- }
- return 0;
- }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
