# Create a 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 <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;

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);

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();

}