# Date Class and 2D Matrix Class

Basics of Object-Oriented Programming

Problem 1:

We gave a few different versions of the DateClass.
a) Modify the setDatefunction to use a switch(m_month) statement that checks if the day is appropriate for the month given. For February assume thatthere are 29 days if the year is divisible by 4 (this is easy to check using the integer mod function(% ) using year%4 which gives the integer remainder of the year divided by 4 so if this is zero,the year is divisible by 4) and otherwise there are 28. In addition add a check to ensure that the year is in the range
of recorded history in the past (going back to ~-3,000) to some anticipated future where this
code will have passed into long obsolescence (say ~3000 to be generous).
To make testing of this function more straightforward, print out warnings rather than using assert
statements (which will stop the program). Also, have your function return a boolean value of
true if the date is successfully set and false if it does not. You should check this returned
value in the main program and print out another warning if the date was not set successfully.
b) Add a constructor to the class that accepts a month, day, and year as arguments but has a
default date of 1 Jan 2000. Your constructor should call your setDatefunction for non-default
values and use an assert statement on setDate’s return value to ensure the date was set
successfully.
c) Change the print() member function into an overload of the <<operator. Note that the
overload of the <<operator will need to be a friend function.

Problem 2:

We created an IntArray class in the notes that illustrated the concept of RAII. The
code for IntArray is posted on OWL in the folder PS3Code, IntArray.cpp.
a) Using IntArray as a template, create a new class DoubArraythat dynamically creates a
double array using the RAII paradigm. In addition to allocating the space, the constructor
should initialize the array elements to zero (You can initialize a dynamically allocated
array with uniform initialization and if you put no numbers in the brackets it is assumed
all entries are zero, e.g. array=new double {}; will dynamically allocate an
array of length 5 with all entries 0).
b) Add a member function L2Norm that finds the Euclidean norm of the array assuming it
is a vector of values (i.e. if this is a vector, return its magnitude, or length). For a long
vector this can be a costly procedure so it is worthwhile to save the result as a member
variable m_normand return this value if the calculation has been previously performed
by an earlier function call. This necessitates another member variable m_normflag
that should be true if the old norm value can be reused and false if it needs to be
recalculated. These member variables should be initialized in the constructor and the
value of m_normflagshould be reset whenever the values of the array are changed.
c) Add an operator overload for the – operator and the assignment operator = that computes
the difference of two arrays as vectors (i.e. if a, b, and c are of the DoubArrayclass
then a statement like c=a-b; should be possible once these two operators have been
d) Further alter your code so that it uses the DoubArrayclass where the length
is dictated by the number of gridpointsN. This should be the only array type used by
your code.Write a function getNto ask the user for this value. Yourimplementation of the implicit Euler’s methods should also be turned into a functionImpEulerthat accepts as arguments N and y where N is the number of grid points andy is a DoubArrayof appropriate length passed by reference and whose memberm_array is set to the initial value 1.0.
e) To test your code we will make use of the functions you created for the DoubArray
class. You need to ensure all parts of your code are tested.

Write another function to fill a DoubArraywith the exact solution (the arguments
should be the same as your ImpEulerfunction). Rather than use the getNfunction for
this test, loop over values of N=10, 20, 30, …100 and in each iteration of the loop find an
exact solution, the implicit Euler solution, the difference between the exact and Euler
solution (using the overloaded – operator), and the Euclidean norm of the difference
(using the L2Norm function). You should output to a file the value of h and the
Euclidean norm of the difference between the exact and Euler approximation and include
a plot of the norm of the difference versus h. Comment on what you find (Does the
difference appear to go to zero as h goes to zero? Is the plot linear?

Problem 3:

This problem will use the RAII concept for matrices. You are reminded that the
code for IntArray is posted on OWL in the folder PS3Code, IntArray.cpp.
a) Using the IntArray class as a template, create a class DoubMatrixthat dynamically
allocates n x m matrices. (the number of rows and number of columns should be
member variables). The matrices should be created as contiguous blocks of memory. The
copy constructor should use assertions to verify that the two matrices have the same size.
b) You have to take arbitrary sized matrices of class DoubMatrix. This should be a friend
function of the DoubMatrixclass that overloads the * operator.
c) You have to take arbitrary sized matrices of class DoubMatrixand vectors of class
DoubArrayfrom part a) of problem 2 above. The vector-matrix and matrix-vector
multiplications should be overloaded * operators that are friend functions of both
DoubMatrixand DoubArray.
As part c) now involves multiple classes, you should separate your file into a main.cpp driver
program that tests various functions, a DoubArray.cpp and DoubMatrix.cpp for the class
implementations and two header files DoubArray.h and DoubMatrix.h.

Inheritance, Virtual Functions, & Templates

Problem 1: In the last assignment, problem 2, you created the DoubArray class and added some
features in problem 2a, 2b, and 2c.

1. a) Add operator overloads to the DoubArray class for the [] operator, providing both a
version for const objects and non-const objects and providing a check on the validity of
the index. The class implementation should be in a separateDoubArray.cpp file and it should havean associated header DoubArray.h.
b) we wrote several implementations of the “DList¬Series ¬Polynomial”
inheritance tree. Alter the Series and Polynomial classes as necessary to inherit from DoubArray
rather than DList You will need to addconstructors to these classes to pass up the size of the array to be initialized in theDoubList constructor.   Problem 3: Use templates to write a single class that can replace our IntArray and DoubArray
classes.

Solution

abstractodesolver.cpp

//problem 2.1

//so writing the methods

#include “abstractodesolver.h”

voidAbstractOdeSolver::SetStepSize(double h)

{

stepSize = h;

}

voidAbstractOdeSolver::SetTimeInterval(double t0, double t1)

{

initialTime = t0;

finalTime = t1;

}

voidAbstractOdeSolver::SetInitialValue(double y0)

{

initialValue = y0;

}

abstractodesolver.h

#ifndef ABSTRACTODESOLVER_H

#define ABSTRACTODESOLVER_H

classAbstractOdeSolver

{

public:

voidSetStepSize(double h);

voidSetTimeInterval(double t0, double t1);

voidSetInitialValue(double y0);

virtual double RightHandSide(double y, double t) = 0;

virtual double SolveEquation() = 0;

protected:

doublestepSize;

doubleinitialTime;

doublefinalTime;

doubleinitialValue;

};

#endif // ABSTRACTODESOLVER_H

doubarray.cpp

#include “doubarray.h”

DoubArray::DoubArray(int length) // constructor

{

assert(length > 0);

m_array = new double[length]{};

m_norm = 0.0; //because 0 init

m_normflag = true; //because 0 init

m_length = length;

//std::cout<< “…finishing constructor\n”;

}

DoubArray::DoubArray(constDoubArray&org_array) // copy constructor

{

m_normflag = org_array.m_normflag;

m_norm = org_array.m_norm;

m_length = org_array.m_length;

//delete first else memory leak !!!!

delete[] m_array;

m_array = new double[m_length];

for (int i=0; i<m_length; i++)

m_array[i]=org_array.m_array[i];

//std::cout<< “…finishing copy constructor\n”;

}

DoubArray::~DoubArray() // destructor

{

// Dynamically delete the array we allocated earlier

delete[] m_array ;

//std::cout<< “…finishing destructor\n”;

}

voidDoubArray::setValue(int index, double value)

{

reset();

m_array[index] = value;

}

voidDoubArray::reset()

{

m_normflag = false;

}

double&DoubArray::getValue(int index)

{

returnm_array[index];

}

intDoubArray::getLength()

{

returnm_length;

}

doubleDoubArray::L2Norm()

{

if(m_normflag)

returnm_norm;

m_norm = 0.0;

for(int i = 0; i <m_length; ++i)

m_norm += std::pow(m_array[i],2.0);

returnm_normflag = true, m_norm;

}

DoubArray&DoubArray::operator =(DoubArrayconst&other)

{

m_normflag = other.m_normflag;

m_norm = other.m_norm;

m_length = other.m_length;

delete[] m_array;

m_array = new double[m_length];

for (int i=0; i<m_length; i++)

m_array[i]=other.m_array[i];

return *this;

}

DoubArrayDoubArray::operator -(DoubArrayconst&other) const

{

DoubArraydoubArray(m_length);

doubArray.reset();

for (int i=0; i<m_length; i++)

doubArray.m_array[i]=m_array[i]-other.m_array[i];

returndoubArray;

}

doubleconst&DoubArray::operator [](int index)const

{

assert((index >= 0) && (m_length> index));

returnm_array[index];

}

double&DoubArray::operator [](int index)

{

assert((index >= 0) && (m_length> index));

reset();

returnm_array[index];

}

doubarray.h

//doubarray class from past assignment

#ifndef DOUBARRAY_H

#define DOUBARRAY_H

#include <iostream>

#include <cmath>

#include <cassert>

classDoubArray

{

private:

intm_length; //array length

doublem_norm; //array discrete L2Norm

boolm_normflag; //array discrete L2Norm validation flag

protected:

double *m_array; //underlying array

public:

DoubArray(int length); //constructor

DoubArray(constDoubArray&org_array); //copy constructor

virtual ~DoubArray(); //virtual destructor

virtual void reset(); //resets flags

voidsetValue(int index, double value); //set array value by index

double&getValue(int index); //get value by index

intgetLength(); //get array length

double L2Norm(); //get l2Norm

DoubArray operator -(DoubArrayconst&other) const; //minus operator overload

doubleconst&operator [](int index)const; //bracket op overload

double&operator [](int index); //bracket op overload

};

#endif // DOUBARRAY_H

forwardeulersolver.cpp

#include “forwardeulersolver.h”

#include <fstream>

voidForwardEulerSolver::SetRightHandSide(rightHandSide_tconst&value)

{

rightHandSide = value;

}

doubleForwardEulerSolver::RightHandSide(double y, double t)

{

returnrightHandSide(y,t);

}

voidForwardEulerSolver::SetFileName(std::string const&value)

{

fileName = value;

}

doubleForwardEulerSolver::SolveEquation()

{

//prepare file to write output

std::ofstreamfos(fileName.c_str());

if(!fos.is_open())

throwstd::runtime_error(“Can’t open file outExpEuler.data”);

//init initial value problem (IVP)

double t(initialTime);

double y(initialValue);

//write IVP to file

fos<< t << ” ” << y << “\n”;

while(t <finalTime)

{

//calculate new stepsize

double h = std::min(stepSize,finalTime-t);

//update solution

y = y + h*RightHandSide(y,t);

//update time node

t = t+h;

//write value snapshot to file

fos<< t << ” ” << y << “\n”;

}

//return final solution

return y;

}

forwardeulersolver.h

//problem 2.2

#ifndef FOWARDEULERSOLVER_H

#define FOWARDEULERSOLVER_H

#include “abstractodesolver.h”

#include <functional>

classForwardEulerSolver :

publicAbstractOdeSolver

{

public:

typedefstd::function<double(double,double)>rightHandSide_t;

doubleRightHandSide(double y, double t);

doubleSolveEquation();

//set right hand side

voidSetRightHandSide(conststd::function<double (double y, double t)>&value);

voidSetFileName(std::string const&value);

private:

rightHandSide_trightHandSide;

std::stringfileName;

};

#endif // FOWARDEULERSOLVER_H

intarray.h

#ifndef INTARRAY_H

#define INTARRAY_H

#include <iostream>

#include <cassert>

classIntArray

{

private:

int *m_array;

intm_length;

public:

IntArray(int length) // constructor

{

assert(length > 0);

m_array = new int[length];

m_length = length;

std::cout<< “…finishing constructor\n”;

}

IntArray(constIntArray&org_array) // copy constructor

{

m_length = org_array.m_length;

m_array = new int[m_length];

for (int i=0; i<m_length; i++)

m_array[i]=org_array.m_array[i];

std::cout<< “…finishing copy constructor\n”;

}

~IntArray() // destructor

{

// Dynamically delete the array we allocated earlier

delete[] m_array ;

std::cout<< “…finishing destructor\n”;

}

voidsetValue(int index, int value) { m_array[index] = value; }

int&getValue(int index) { return m_array[index]; }

intgetLength() { return m_length; }

};

#endif // INTARRAY_H

laplacesolver.h

//able to calculate potential (EvalSum) given radius r and angle theta

#ifndef LAPLACESOLVER_H

#define LAPLACESOLVER_H

#include “rthetalegendre.h”

classLaplaceSolver

{

public:

typedef double(*V_t)(double);

private:

double *angles; //array containing theta’s

doubledTheta; //stepsize in theta

int M; //angles size

V_t V; //pased user defined function

public:

//constructor

LaplaceSolver(int N,

int M,

V_t V)

{

//allocmemmyRTheta

myRTheta = new RThetaLegendre(N);

//discrete version of coefficients is not correct in task 1e)

//integral should be approxed from 0 to pi with M points

//so this gives stepsize h[i] = i*pi/(M-1) so h = 0 and h[M-1] = pi

//1) make size M array

angles = new double[M];

//stepsizedTheta is pi/(size(angles)-1) !!!

dTheta = M_PI/(double(M-1));

//populate angles correct formula theta = i*pi/(M-1) = i*dTheta

for(int i = 0; i < M; ++i)

angles[i] = double(i)*dTheta;

this->M = M;

this->V = V;

}

~LaplaceSolver()

{

deletemyRTheta;

delete[] angles;

}

//evaluate sum

doubleEvalSum(double r, double theta)

{

returnmyRTheta->EvalSum(r,theta);

}

//evaluates the coefficients

voidEvalCoefficients()

{

//bluff because r^(-n-1)  = 1.0 if r = 1.0

myRTheta->SetR(1.0);

for(int i = 0; i < M; ++i) //start with integral loop in order to keep x valid in polynomial calc

{

double theta = angles[i];

doublesinTheta = std::sin(theta);

doublecosTheta = std::cos(theta);

doubleVTheta = V(theta);

doubletmp = VTheta*sinTheta*dTheta;

for(int n = 0; n <myRTheta->getLength(); ++n)

(*myRTheta)[n] += tmp*(n+1.0/2.0)*myRTheta->EvalPoly(n,cosTheta);

}

}

};

#endif // LAPLACESOLVER_H

legendre.h

#ifndef LEGENDRE_H

#define LEGENDRE_H

#include “series.h”

class Legendre :

public Series

{

public:

doublem_lft;

doublem_rgt;

doublem_h;

Legendre(int N,

doublem_lft = -1.0,

doublem_rgt = +1.0) :

Series(N)

{

this->m_lft = m_lft;

this->m_rgt = m_rgt;

m_h = (this->m_rgt – this->m_lft)/double(N-1);

}

virtual ~Legendre()

{

}

//evaluate polynomial (warning x must stay same until reset triggered with n = 0)

virtual double EvalPoly(int n, double x)

{

if (!((x <= m_rgt) && (x >=m_lft)))

std::cout<< “Warning: Outside expected range for this polynomial. \n”;

static double Pnm1, Pn, myX;

if(n == 0)

{

myX = x;

returnPn = 1.0;

}

else if(n > 0)

{

if(myX != x)

{

std::cout<< “warning x must stay same until reset triggered with n = 0”;

exit(1);

}

//substite n -> n-1 from sheet formula

Pnm1 = ((2.0*n-1.0)*x*Pn – (n-1.0)*Pnm1)/double(n);

//swap pnm1, pn

std::swap(Pnm1,Pn);

returnPn;

}

//not allowed to reach if so return nan

return 1.0/0.0;

}

doubleEvalSum(double x)

{

doubletmp(0.0);

for(int i = 0; i <getLength(); ++i)

tmp+=m_array[i]*EvalPoly(i,x);

returntmp;

}

};

#endif // LEGENDRE_H

main.cpp

#include <iostream>

#include <fstream>

#include <cmath>

#include “intarray.h”

#include “doubarray.h”

#include “marray.h”

#include “polynomial.h”

#include “legendre.h”

#include “rthetalegendre.h”

#include “laplacesolver.h”

#include <sstream>

#include “forwardeulersolver.h”

#include “rungekuttasolver.h”

//problem 1a)

voiddoubArrayBracketOperatorDemo()

{

//test by simply check if c = a-b

DoubArraya(1), b(1), c(1);

a = 5.0;

b = 3.0;

c = a-b;

assert(c == 2.0);

std::cout<< “Problem 1a test passed” <<std::endl;

}

//problem 1b)

voidinheritDoubArrayDemo()

{

//analogous test as for problem 1a)

{

Polynomial a(1), b(1), c(1);

a = 5.0;

b = 3.0;

c = a-b;

assert(c == 2.0);

}

//analogous test as for problem 1a)

{

Series a(1), b(1), c(1);

a = 5.0;

b = 3.0;

c = a-b;

assert(c == 2.0);

}

std::cout<< “Problem 1b test passed” <<std::endl;

}

//problem 1c)

voidlegendrePolyDemo()

{

//open for output

std::ofstreamfos(“legendre5.dat”);

//make Legendre series with 30 polynomials

Legendre legendre(30);

//evaluate polynomials

for(double x  = legendre.m_lft; x <legendre.m_rgt; x+= legendre.m_h)

{

fos<< x << ” “;

for(int i = 0; i < 5; ++i)

{

legendre[i] = 1.0;

fos<<legendre.EvalSum(x) << ” “;

legendre[i] = 0.0;

}

fos<< “\n”;

}

std::cout<< “Problem 1c legendre5.dat written (plot by e.g. gnuplot plot.gp)” <<std::endl;

}

//problem 1d)

voidrThetaLegendreDemo()

//an = 1/n not possible due to singularity @ n = 0

//so i chose an = 1/n, and a0 = 1

{

//open for output

std::ofstreamfos(“rThetaLegendre4.dat”);

//set an = 1/(n+1)

RThetaLegendrerThetaLegendre(100);

rThetaLegendre = 1.0;

for(int i = 1; i <rThetaLegendre.getLength(); ++i)

rThetaLegendre[i] = 1.0/double(i);

//evaluate series

for(double theta  = 0.0; theta < M_PI; theta+= M_PI/99.0)

{

fos<< theta << ” “;

fos<<rThetaLegendre.EvalSum(1,theta) << ” “;

fos<<rThetaLegendre.EvalSum(2,theta) << ” “;

fos<<rThetaLegendre.EvalSum(4,theta) << ” “;

fos<<rThetaLegendre.EvalSum(8,theta) << ” “;

fos<< “\n”;

}

std::cout<< “Problem 1d rThetaLegendre4.dat written (plot by e.g. gnuplot plot.gp)” <<std::endl;

}

//problem 1e)

voidlaplaceSolverDemo()

{

//open for output

std::ofstreamfos(“laplaceSolver4.dat”);

LaplaceSolverlaplaceSolver(/*evenly spaced points [0,pi]*/100,

/*numPoints for integral approximation*/1e4,

/*V(theta)=*/[](double theta){ return std::sin(2.0*theta);});

//set coefficients

laplaceSolver.EvalCoefficients();

//evaluate series

for(double theta  = 0.0; theta < M_PI; theta+= M_PI/99.0)

{

fos<< theta << ” “;

fos<<laplaceSolver.EvalSum(1,theta) << ” “;

fos<<laplaceSolver.EvalSum(2,theta) << ” “;

fos<<laplaceSolver.EvalSum(4,theta) << ” “;

fos<<laplaceSolver.EvalSum(8,theta) << ” “;

fos<< “\n”;

}

std::cout<< “Problem 1e laplaceSolver4.dat written (plot by e.g. gnuplot plot.gp)” <<std::endl;

}

//problem 2)

{

//lambda the RightHandSide given

autoRightHandSide = [](double /*y*/, double t) { return 1.0 + t;};

//lambda the exact solution given

autoExactSolution = [](double t){ return 0.5*(t*t+2*t+4.0);};

//init explicit euler

ForwardEulerSolverforwardEulerSolver;

RungeKuttaSolverrungeKuttaSolver;

//initial values for IVP

doubleinitialTime(0.0);

doublefinalTime(1.0);

doubleinitialValue(2.0);

doublestepSize(0.1);

//init solver

forwardEulerSolver.SetInitialValue(initialValue);

forwardEulerSolver.SetRightHandSide(RightHandSide);

forwardEulerSolver.SetTimeInterval(initialTime,finalTime);

forwardEulerSolver.SetStepSize(stepSize);

rungeKuttaSolver.SetInitialValue(initialValue);

rungeKuttaSolver.SetRightHandSide(RightHandSide);

rungeKuttaSolver.SetTimeInterval(initialTime,finalTime);

rungeKuttaSolver.SetStepSize(stepSize);

//1) check solution, by checking relative error

//   ||yApprox – yExact||/||yExact|| <rtol

{

//set output file name

forwardEulerSolver.SetFileName(“ForwardEulerSolver_h01Run.dat”);

rungeKuttaSolver.SetFileName(“RungeKuttaSolver_h01Run.dat”);

double yApprox1  = forwardEulerSolver.SolveEquation();

std::cout<< “ForwardEulerSolver_h01Run.dat written” <<std::endl;

double yApprox2  = rungeKuttaSolver.SolveEquation();

std::cout<< “RungeKuttaSolver_h01Run.dat written” <<std::endl;

double relErr1 = std::abs(yApprox1 – ExactSolution(1.0))/std::abs(ExactSolution(1.0));

double relErr2 = std::abs(yApprox2 – ExactSolution(1.0))/std::abs(ExactSolution(1.0));

std::cout<< “ForwardEulerSolver with stepsize 0.001 has relative error = ” << relErr1 <<std::endl;

std::cout<< “RungeKuttaSolver with stepsize 0.001 has relative error = ” << relErr2 <<std::endl;

}

//2) how the choice of the size affects the accuracy ?

// error goes like E ~ h ^ p, with E:error, h:stepSize and p order of method

// so E1/E2 = (h1/h2)^p ~> p = log(E1/E2)/log(h1/h2)

// now expecting p = 1 for euler and p = 4 for rk4

// but rightHandSide only scales with t so we have exact solution with all integrators where p > 1

// therefore for rk4 we have independant of stepSize the exact solution in this case (apart from possible instabilities)

{

//estimate order for euler method

double h1 = 0.5*(finalTime-initialTime);

forwardEulerSolver.SetStepSize(h1);

forwardEulerSolver.SetFileName(“ForwardEulerSolver_h05Run.dat”);

double yApprox1  = forwardEulerSolver.SolveEquation();

std::cout<< “ForwardEulerSolver_h05Run.dat written” <<std::endl;

double h2 = (finalTime-initialTime);

forwardEulerSolver.SetStepSize(h2);

forwardEulerSolver.SetFileName(“ForwardEulerSolver_h10Run.dat”);

double yApprox2 = forwardEulerSolver.SolveEquation();

std::cout<< “ForwardEulerSolver_h10Run.dat written” <<std::endl;

double p = std::log(std::abs(yApprox1-ExactSolution(finalTime))/std::abs(yApprox2-ExactSolution(finalTime)))/std::log(h1/h2);

std::cout<< “Estimated ForwardEulerSolver order is “<< p <<std::endl;

}

}

//problem 3)

voidmArrayDemo()

{

//test by comparing sum of Array<int> Array<double>

int n = 50;

Array<int>iArray(n);

Array<double>dArray(n);

//initi

for(int i = 0; i < n; ++i)

{

iArray[i] = 4;

dArray[i] = 4.5;

}

//calculate sum

intiSum(0);

doubledSum(0);

for(int i = 0; i < n; ++i)

{

iSum += iArray[i];

dSum += dArray[i];

}

//check if correct

assert(iSum == 50*4);

assert(std::abs(dSum – 4.5*50.0) < 1.e-6);

std::cout<< “Problem 3 test passed” <<std::endl;

}

int main(int , char *[])

{

std::cout<< “The beginning” <<std::endl;

std::cout<< “\nProblem1” <<std::endl;

doubArrayBracketOperatorDemo();

inheritDoubArrayDemo();

legendrePolyDemo();

rThetaLegendreDemo();

laplaceSolverDemo();

//problem2)

std::cout<< “\nProblem2” <<std::endl;

std::cout<< “\nProblem3” <<std::endl;

mArrayDemo();

std::cout<< “\nThe end” <<std::endl;

return 0;

}

marray.h

//Problem 3)

//basically replaced all double inside DoubArray with template argument _T

//keep in mind nobody would use these arrays in real life…

//nowadays expression templates are used in state of the art implementations

#ifndef MARRAY_H

#define MARRAY_H

template<typename _T>

class Array

{

private:

intm_length;

_T m_norm;

boolm_normflag;

protected:

_T *m_array;

public:

Array(int length) // constructor

{

assert(length > 0);

m_array = new _T[length]{};

m_norm = 0.0; //because 0 init

m_normflag = true; //because 0 init

m_length = length;

//std::cout<< “…finishing constructor\n”;

}

Array(const Array &org_array) // copy constructor

{

m_normflag = org_array.m_normflag;

m_norm = org_array.m_norm;

m_length = org_array.m_length;

//delete first else memory leak !!!!

delete[] m_array;

m_array = new _T[m_length];

for (int i=0; i<m_length; i++)

m_array[i]=org_array.m_array[i];

//std::cout<< “…finishing copy constructor\n”;

}

~Array() // destructor

{

// Dynamically delete the array we allocated earlier

delete[] m_array ;

//std::cout<< “…finishing destructor\n”;

}

virtual void reset()

{

m_normflag = false;

}

voidsetValue(int index, _T value)

{

reset();

m_array[index] = value;

}

_T&getValue(int index)

{

returnm_array[index];

}

intgetLength()

{

returnm_length;

}

_T L2Norm()

{

if(m_normflag)

returnm_norm;

m_norm = 0.0;

for(int i = 0; i <m_length; ++i)

m_norm += std::pow(m_array[i],2.0);

returnm_normflag = true, m_norm;

}

Array &operator =(Array const&other)

{

m_normflag = other.m_normflag;

m_norm = other.m_norm;

m_length = other.m_length;

delete[] m_array;

m_array = new _T[m_length];

for (int i=0; i<m_length; i++)

m_array[i]=other.m_array[i];

return *this;

}

Array operator -(Array const&other)

{

Array<_T>mArray(m_length);

mArray.reset();

for (int i=0; i<m_length; i++)

mArray.m_array[i]=m_array[i]-other.m_array[i];

returnmArray;

}

_T const&operator [](int index)const

{

assert((index >= 0) && (m_length> index));

returnm_array[index];

}

_T &operator [](int index)

{

assert((index >= 0) && (m_length> index));

reset();

returnm_array[index];

}

};

#endif // MARRAY_H

polynomial.h

#ifndef POLYNOMIAL_H

#define POLYNOMIAL_H

#include <iostream>

#include <cassert>

#include <cmath>

#include “series.h”

class Polynomial :

public Series

{

private :

doublem_lft=-1.;

doublem_rgt=1.;

public:

Polynomial(double leftendpoint=-1.0,

doublerightendpoint=1.0,

int N=10) : Series(N)

{

m_lft=leftendpoint;

m_rgt=rightendpoint;

//std::cout<< “…constructing a Polynomial” <<std::endl;

}

doubleEvalPoly(double x)

{

if (!((x <m_rgt) && (x>=m_lft)))

std::cout<< x << “Warning: Outside range for this polynomial. \n”;

doublepx=0., xn=1.0;

for (int n = 0; n <getLength(); n++)

{

px += m_array[n]*xn;

xn *= x;

}

returnpx;

}

virtual ~Polynomial()

{

}

bool decreasing(const double x)

{

intmaxN=getLength();

doublexn=1.0;

for (int i = 1; i <maxN; i++)

if (fabs(m_array[i-1]*xn) <fabs(m_array[i]*(xn*x)))

{

xn=xn*x;

return false;

}

return true;

}

};

#endif // POLYNOMIAL_H

rthetalegendre.h

#ifndef RTHETALEGENDRE_H

#define RTHETALEGENDRE_H

#include “legendre.h”

classRThetaLegendre :

//legendre base

public Legendre

{

private:

//use hint and make r member

double r;

public:

//use base constructor

using Legendre::Legendre;

//make virtual destructor

virtual ~RThetaLegendre() { }

//Evaluate polynoms

doubleEvalPoly(int n, double cosTheta)

//polynoms are pow(r,-(n+1))*Pn(cos(theta))

{

returnstd::pow(r,-(n+1.0))*Legendre::EvalPoly(n,cosTheta);

}

//set r

voidSetR(double value)

{

this->r = value;

}

//evaluate sum

doubleEvalSum(double r, double theta)

{

//set r

SetR(r);

//use base’s evalsum

return Legendre::EvalSum(std::cos(theta));

}

};

#endif // RTHETALEGENDRE_H

rungekuttasolver.h

//problem 2

#ifndef RUNGEKUTTASOLVER_H

#define RUNGEKUTTASOLVER_H

#include “abstractodesolver.h”

#include <functional>

classRungeKuttaSolver :

publicAbstractOdeSolver

{

public:

typedefstd::function<double(double,double)>rightHandSide_t;

doubleRightHandSide(double y, double t);

doubleSolveEquation();

voidSetRightHandSide(conststd::function<double (double y, double t)>&value);

voidSetFileName(std::string const&value);

private:

rightHandSide_trightHandSide;

std::stringfileName;

};

#endif // RUNGEKUTTASOLVER_H

series.h

#ifndef SERIES_H

#define SERIES_H

#include <iostream>

#include <cassert>

#include <cmath>

#include “doubarray.h”

class Series :

publicDoubArray

{

private :

doublem_sum;

boolm_sumflag;

public:

Series(int length) :

DoubArray(length)

{

m_sum = 0.0;

reset();

}

virtual ~Series()

{

}

virtual void reset()

{

DoubArray::reset(); //first reset base

m_sumflag = false;

}

bool decreasing()

{

intmaxN=getLength();

for (int i = 1; i <maxN; i++)

if (fabs(m_array[i]) >= fabs(m_array[i-1]))

return false;

return true;

}

doublegetSum()

{

if (m_sumflag)

returnm_sum;

m_sum=0.0;

intmaxN=getLength();

for (int i =0; i <maxN; i++)

m_sum += m_array[i];

m_sumflag = true;

returnm_sum;

}

};

#endif // SERIES_H