Instructions
Requirements and Specifications
Source Code
#include <iostream>
#include <cmath>
#include <string>
#include <sstream>
#include <fstream>
#include <iomanip>
#include <vector>
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<V0> 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<V0> 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<V0> 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<V0> forward_difference(float Vinit, float Vs, float dt)
{
// First, create the vector
vector<V0> 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<V0> backward_difference(float Vinit, float Vs, float dt)
{
// First, create the vector
vector<V0> 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<V0> central_difference(float Vinit, float Vs, float dt)
{
// First, create the vector
vector<V0> 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<V0> V0sol_forward = forward_difference(Vinit, Vs, dt);
vector<V0> V0sol_backward = backward_difference(Vinit, Vs, dt);
vector<V0> 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"<<endl;
// 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();
}