Program to Solve ODE Equation in C++ Assignment Solution

Instructions

Objective

Write a program to solve ODE equation in C++ language.

Requirements and Specifications

Source Code

#include #include #include #include #include #include #include using namespace std; // define the size of arrays here. This is also the number of points #define SIZE 5000 // We will create a class that will be useful to store the solution at each time step class V0 { private: float time; // time value float value; // value of voltage public: V0() { time = 0; value = 0; }; V0(float v, float t) { time = t; value = v; } // getters float getTime() {return time;} float getValue() {return value;} }; V0 Vnew_forward(vector V0arr, float Vs, float dt, int i) { /* Given the index i, the vector V and time-step dt, this function calculates the new Value at the inex i+1 according to the Forward Difference method This function basically calculates the value of the Voltage at index i+1 (or x +dx) */ float input = 5*Vs; float V0i_term = V0arr[i].getValue()*(5 - 2/pow(dt,2) - 2/dt); // i term float V0im1_term = V0arr[i-1].getValue()*(1/pow(dt,2)); // i-1 term float divisor = 1/pow(dt,2) + 2/dt; float val = (input - V0i_term - V0im1_term)/divisor; return V0(val, (i+1)*dt); } V0 Vnew_backward(vector V0arr, float Vs, float dt, int i) { /* Given the index i, the vector V and time-step dt, this function calculates the new Value at the inex i+1 according to the Central Difference method This function basically calculates the value of the Voltage at index i+1 (or x +dx) */ float input = 5*Vs; float V0i_term = V0arr[i].getValue()*(5 - 2/pow(dt,2) + 2/dt); // i term float V0im1_term = V0arr[i-1].getValue()*(1/pow(dt,2) - 2/dt); // i-1 term float val = (input - V0i_term - V0im1_term)*pow(dt,2); return V0(val, (i+1)*dt); } V0 Vnew_central(vector V0arr, float Vs, float dt, int i) { /* Given the index i, the vector V and time-step dt, this function calculates the new Value at the inex i+1 according to the Forward Difference method This function basically calculates the value of the Voltage at index i+1 (or x +dx) */ float input = 5*Vs; float V0i_term = V0arr[i].getValue()*(5 - 2/pow(dt,2)); // i term float V0im1_term = V0arr[i-1].getValue()*(1/pow(dt,2) - 1/dt); // i-1 term float divisor = 1/pow(dt,2) + 1/dt; float val = (input - V0i_term - V0im1_term)/divisor; return V0(val, (i+1)*dt); } vector forward_difference(float Vinit, float Vs, float dt) { // First, create the vector vector Vsol; Vsol.push_back(V0(Vinit, 0)); Vsol.push_back(V0(Vinit, dt)); for(int i = 1; i < SIZE-1; i++) { V0 V0new = Vnew_forward(Vsol, Vs, dt, i); Vsol.push_back(V0new); } return Vsol; } vector backward_difference(float Vinit, float Vs, float dt) { // First, create the vector vector Vsol; Vsol.push_back(V0(Vinit, 0)); Vsol.push_back(V0(Vinit, dt)); for(int i = 1; i < SIZE-1; i++) { V0 V0new = Vnew_backward(Vsol, Vs, dt, i); Vsol.push_back(V0new); } return Vsol; } vector central_difference(float Vinit, float Vs, float dt) { // First, create the vector vector Vsol; Vsol.push_back(V0(Vinit, 0)); Vsol.push_back(V0(Vinit, dt)); for(int i = 1; i < SIZE-1; i++) { V0 V0new = Vnew_central(Vsol, Vs, dt, i); Vsol.push_back(V0new); } return Vsol; } float Vanalytic(float t) { /* Return the analytic solution at time t */ float V =10 - exp(-t)*(10*cos(2*t) + 5*sin(2*t)); return V; } int main() { // Variable to store the name of the output file string file_name; // Ask for filename cout << "Enter output file: "; getline(cin, file_name); // Define parameters here like input value (source), initial condition, initial time, final time and time step float Vinit = 0; // V at t = 0 float V0p = 0; // first derivative of V at t = 0 float Vs = 10; // Define the final time of simulation float tf = 5.0; // in seconds // Calculate step float dt = tf/(float)SIZE; // Solve using the different methods vector V0sol_forward = forward_difference(Vinit, Vs, dt); vector V0sol_backward = backward_difference(Vinit, Vs, dt); vector V0sol_central = central_difference(Vinit, Vs, dt); // Create array to store time float *Van = new float[SIZE]; float *t= new float[SIZE]; // Calculate all the values for the next time-steps t[0] = 0; for(int i = 1; i < SIZE; i++) { // Store also the time values t[i] = t[i-1] + dt; } // open output file ofstream oFile(file_name); // Print headers cout << left << setw(10) << "Time" << right << setw(20) << "Backward" << setw(20) << "Forward" << setw(20) << "Central" << setw(20) << "V0 - analytical" << endl; cout << left << setw(10) << "----" << right << setw(20) << "--------" << setw(20) << "-------" << setw(20) << "-------" << setw(20) << "---------------" << endl; // Write header to output file //oFile << "time,Vsolution,Vanalytical"< // Loop through all the values and append to the file // Print every 100 values to the console // Create variables to store maximum and minimum voltages, and their times V0 V0max_forward; V0 V0max_backward; V0 V0max_central; V0 V0min_forward; V0 V0min_backward; V0 V0min_central; // Define here the number of values at which the lines will be printed int Nlines = 100; for(int i = 0; i < SIZE; i++) { Van[i] = Vanalytic(t[i]); // Calculate the analytic solution for time i // Write output file as a csv file oFile << t[i] << "," << V0sol_forward[i].getValue() << "," << V0sol_backward[i].getValue() << "," << V0sol_central[i].getValue() << "," << Van[i] << endl; if(i%Nlines == 0) {// Print only every 100 values cout << left << setw(10) << t[i] << right << setw(20) << V0sol_forward[i].getValue() << setw(20) << V0sol_backward[i].getValue() << setw(20) << V0sol_central[i].getValue() << setw(20) << Van[i] << endl; } if(V0sol_forward[i].getValue() < V0min_forward.getValue()) // voltage i is lower than the minimum voltage registered { V0min_forward = V0sol_forward[i]; } if(V0sol_backward[i].getValue() < V0min_backward.getValue()) // voltage i is lower than the minimum voltage registered { V0min_backward = V0sol_backward[i]; } if(V0sol_central[i].getValue() < V0min_central.getValue()) // voltage i is lower than the minimum voltage registered { V0min_central = V0sol_central[i]; } if(V0sol_forward[i].getValue() > V0max_forward.getValue()) // voltage i is lower than the minimum voltage registered { V0max_forward = V0sol_forward[i]; } if(V0sol_backward[i].getValue() > V0max_backward.getValue()) // voltage i is lower than the minimum voltage registered { V0max_backward = V0sol_backward[i]; } if(V0sol_central[i].getValue() > V0max_central.getValue()) // voltage i is lower than the minimum voltage registered { V0max_central = V0sol_central[i]; } } // Display max and min cout << endl << "The maximum voltage obtained using Forward Difference Method is " << V0max_forward.getValue() << " (V) and it happened at t = " << V0max_forward.getTime() << " (s)" << endl; cout << "The minimum voltage obtained using Forward Difference Method is " << V0min_forward.getValue() << " (V) and it happened at t = " << V0min_forward.getTime() << " (s)" << endl; cout << endl << "The maximum voltage obtained using Backward Difference Method is " << V0max_backward.getValue() << " (V) and it happened at t = " << V0max_backward.getTime() << " (s)" << endl; cout << "The minimum voltage obtained using Backward Difference Method is " << V0min_backward.getValue() << " (V) and it happened at t = " << V0min_backward.getTime() << " (s)" << endl; cout << endl << "The maximum voltage obtained using Central Difference Method is " << V0max_central.getValue() << " (V) and it happened at t = " << V0max_central.getTime() << " (s)" << endl; cout << "The minimum voltage obtained using Central Difference Method is " << V0min_central.getValue() << " (V) and it happened at t = " << V0min_central.getTime() << " (s)" << endl; oFile.close(); // Finally delete the vectors delete[] t; delete[] Van; V0sol_backward.clear(); V0sol_forward.clear(); V0sol_central.clear(); }