Commit 1db6be6b authored by Michele Invernizzi's avatar Michele Invernizzi
Browse files

added solutions ex11

parent 7f290b9d
cmake_minimum_required(VERSION 2.8)
project(SimpsonIntegrationUsingAbstractClass)
set(CMAKE_CXX_FLAGS "-Wall -std=c++11")
# Create executable
add_executable(factory factory.cpp)
# Install executable
set(CMAKE_INSTALL_PREFIX ${CMAKE_HOME_DIRECTORY}/install)
install(TARGETS factory DESTINATION bin)
/* Programming Techniques for Scientific Simulations, HS 2016
* Exercise 11.2
*/
#include <iostream>
#include <cstdlib>
#include <cmath>
#include <string>
#include <stdexcept>
// abstract base class version using inheritance and runtime polymorphism
struct Function {
typedef double result_type;
typedef double argument_type;
virtual result_type operator()(argument_type x) const=0;
virtual ~Function() {}
};
// Simpson integration routine
Function::result_type integrate(const Function::argument_type a, const Function::argument_type b, const unsigned bins, const Function& f) {
if (bins==0)
throw std::logic_error("integrate(..) : number of bins has to be positive.");
typedef Function::argument_type arg_type;
typedef Function::result_type res_type;
arg_type half_bin_size=(b-a)/static_cast<arg_type>(2*bins), x_;
res_type sum_odd=f(a+half_bin_size), sum_even(0);
for (unsigned i=1; i<bins; ++i) {
x_=a+2*i*half_bin_size;
sum_even+=f(x_);
sum_odd+=f(x_+half_bin_size);
}
return (2.*sum_even + 4.*sum_odd + f(a) + f(b)) * half_bin_size / 3.;
}
struct sin_lambda : public Function
{
sin_lambda( argument_type lambda = 1. ): lambda_(lambda){}
result_type operator()( argument_type x ) const { return std::sin(lambda_*x); }
private:
argument_type lambda_;
};
struct cos_lambda : public Function
{
cos_lambda( argument_type lambda = 1. ): lambda_(lambda){}
result_type operator()( argument_type x ) const { return std::cos(lambda_*x); }
private:
argument_type lambda_;
};
struct exp_lambda : public Function
{
exp_lambda( argument_type lambda = 1. ): lambda_(lambda){}
result_type operator()( argument_type x ) const { return std::exp(lambda_*x); }
private:
argument_type lambda_;
};
Function* integrandFactory(std::string algorithm , Function::argument_type lambda )
{
if (algorithm=="sin_lambda") {
return new sin_lambda( lambda );
} else if (algorithm=="cos_lambda") {
return new cos_lambda( lambda );
} else if (algorithm=="exp_lambda") {
return new exp_lambda( lambda );
} else {
throw std::invalid_argument("Integrand must be one of the following: 'sin_lambda', 'cos_lambda', 'exp_lambda'.");
}
}
int main(int argc, char ** argv)
{
if (argc != 5 && argc != 6 ) {
std::cerr << argv[0] << " <a> <b> <bins> <function_name> <lambda (optional)> " << std::endl;
exit(1);
}
double a = std::stod( argv[1] );
double b = std::stod( argv[2] );
unsigned bins = std::stoi( argv[3] );
double lambda = 1.;
if (argc == 6)
lambda = std::stod( argv[5] );
Function* myintegrand = integrandFactory( std::string( argv[4] ) , lambda );
std::cout.precision(15);
std::cout << "Integral is: " << integrate( a , b , bins , *myintegrand ) << std::endl;
delete myintegrand;
return 0;
}
cmake_minimum_required(VERSION 2.8)
# you may manually set the compiler if you want
#SET(CMAKE_CXX_COMPILER "c++")
# set the flags for the compiler: all warnings
SET(CMAKE_CXX_FLAGS "-Wall -std=c++11")
SET(CMAKE_BUILD_TYPE RELEASE)
if(APPLE)
find_library(ACCELERATE_LIBRARY Accelerate )
mark_as_advanced(ACCELERATE_LIBRARY)
set(EXTRA_LIBS ${ACCELERATE_LIBRARY})
else()
set(EXTRA_LIBS lapack blas gfortran)
endif()
set(PROGRAM harmonic_chain)
add_executable(${PROGRAM} ${PROGRAM}.cpp)
target_link_libraries(${PROGRAM} ${EXTRA_LIBS})
set(CMAKE_INSTALL_PREFIX ${CMAKE_HOME_DIRECTORY}/install)
install(TARGETS ${PROGRAM} DESTINATION bin)
/* Programming Techniques for Scientific Simulations, HS 2016
* Exercise 11.0
*/
#include <iostream>
#include <cmath>
#include <algorithm>
#include <vector>
#include <fstream>
#define _USE_MATH_DEFINES //for M_PI
typedef double calc_type;
typedef int size_type;
extern "C" {
void dsyev_( const char* JOBZ , const char* UPLO, const int* N, double* A, const int* LDA, double* W,
double* WORK, const int* LWORK, int* INFO ) ;
}
inline void print(const std::vector<calc_type>& A, const size_type n,size_type n_show) {
if (n_show>n)
n_show=n;
for (int i=0; i<n_show; ++i) {
for (int j=0; j<n_show; ++j) {
std::cout << A[i*n+j] << " ";
}
if (n>n_show)
std::cout<<"...";
std::cout << std::endl;
}
if (n>n_show)
std::cout << "..." << std::endl;
}
inline void print(const std::vector<calc_type>& A, const size_type n) { print(A,n,n); }
int main() {
// parameters of the task
const calc_type K=1.;
const calc_type mass=1.;
const size_type N=16;
const size_type size_of_workspace=(N+2)*N+1;
size_type N_show=8;
if (N_show>N)
N_show=N;
// data structures
std::vector<calc_type> M(N*N,0);
std::vector<calc_type> omega(N);
std::vector<calc_type> exact_omega(N);
std::vector<calc_type> workspace(size_of_workspace);
// fill in the Hamiltonian
for (size_type n=0; n<N; ++n) { // diagonal
M[n*N+n]=2;
}
for (size_type n=0; n<N-1; ++n) { // diagonal shifted 1 rows above
M[n*N+n+1]=-1;
}
//print to screen a preview
std::cout<<"\nHalf Hamiltonian"<<std::endl;
print(M,N,N_show);
std::cout<<std::endl;
//calculate
size_type info;
const char task='V', triangle='L';
dsyev_(&task, &triangle, &N, &(M.front()), &N, &(omega.front()), &(workspace.front()), &size_of_workspace, &info);
if (info!=0) {
std::cerr<<"An error occured in the dsyev_ function. INFO="<<info<<std::endl;
return -1;
}
const calc_type sqrt_Km=std::sqrt(K/mass);
for (size_type n=0; n<N; ++n)
omega[n]=sqrt_Km*std::sqrt(omega[n]);
//calculate the exact eigenfrequencies
for (size_type n=0; n<N; ++n)
exact_omega[n]=sqrt_Km*std::sqrt(2-2*std::cos(M_PI*(n+1)/(N+1)));
//print a preview to screen
std::cout<<"Omega Eigenvector"<<std::endl;
for (size_type n=0; n<N_show; ++n) {
std::cout<<omega[n]<<" [";
for (size_type i=0; i<N_show; ++i) {
std::cout<<M[n*N+i]<<(i==N-1?"":", ");
}
if (N>N_show)
std::cout<<"...";
std::cout<<"]"<<std::endl;
}
if (N>N_show)
std::cout<<"... [...]"<<std::endl;
std::cout<<std::endl;
//print to file full results
std::ofstream outfile;
outfile.open("results_cpp.data");
outfile<<"# exact_omega omega eigenvector"<<std::endl;
outfile.precision(15);
for (size_type n=0; n<N; ++n)
{
outfile<<exact_omega[n]<<" "<<omega[n];
for (size_type i=0; i<N; ++i)
outfile<<" "<<M[n*N+i];
outfile<<std::endl;
}
return 0;
}
#!/usr/bin/env python
# -*- coding: utf-8 -*-
#
# Programming Techniques for Scientific Simulations, HS 2016
# Exercise 11.1
# File: harmonic_chain.py
from __future__ import (absolute_import, division, print_function,
unicode_literals)
import json
import numpy as np
def main(K=1., m=1., N=16):
H = 2 * np.eye(N) - np.eye(N, k=1)
eigval, eigvec = np.linalg.eigh(H, UPLO='U')
omega = np.sqrt(K / m * eigval)
exact_omega = np.sqrt(
K / m * (2 - 2 * np.cos(
np.pi * np.arange(1, N + 1) / (N + 1)
))
)
# print a preview of the results to screen
eigvec = np.transpose(eigvec) # we want to print eigenvectors on rows
np.set_printoptions(threshold=100)
print(
"\nHalf Hamiltonian:", H,
"\nEigenfrequencies omega:", omega,
"\nEigenvectors:", eigvec,
sep='\n',
end='\n\n'
)
# print the full results to file
np.savetxt(
"results_py.data",
np.c_[exact_omega, omega, eigvec],
header="exact_omega omega eigenvector",
fmt=str('%.15g')
)
# # print to file in a nicer format
# with open('results_py.json', 'w') as f:
# json.dump(
# {
# 'exact_omega': exact_omega.tolist(),
# 'omega': omega.tolist(),
# 'eigvec': eigvec.tolist()
# },
# f
# )
if __name__ == "__main__":
main()
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment