Commit 0f473bee authored by rworreby's avatar rworreby

Add ex 13 solution

parent d63a8085
p=print
a=b=''
i=26
while i:l=chr(123-i);a+=l;b+=l.upper()+l;i-=1
c='-'*23
p(b[::2]+a);p(b);p(c)
while b:p('| %-20s|'%' '.join(b[:10]));b=b[10:]
p(c)
d='prelpamilpoveene'
for e in a:
if e in'aglo':x=e+d[i::4];b+=', %s begins with '%x+e;e=e+': ',x;i+=1
p([*e])
p(b[2:])
\ No newline at end of file
#!/usr/bin/env python
# -*- coding: utf-8 -*-
#
# Programming Techniques for Scientific Simulations, HS 2020
# Exercise 13.1
from __future__ import (absolute_import, division, print_function,
unicode_literals)
import h5py
import numpy as np
def main(K=1.0, m=1.0, N=4):
# Set up H (only upper triangle).
H = 2 * np.eye(N) - np.eye(N, k=1)
# Compute eigenvalues and eigenvectors.
eigval, eigvec = np.linalg.eigh(H, UPLO='U')
print("EW ", eigval)
# Compute frequencies from eigenvalues.
omega = np.sqrt(K / m * eigval)
# Exact frequencies.
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(
"\nHamiltonian (upper triangle):", H,
"\nEigenfrequencies omega:", omega,
"\nEigenvectors:", eigvec,
sep='\n',
end='\n\n'
)
# Print the full results to a file.
np.savetxt(
"results_py.data",
np.c_[exact_omega, omega, eigvec],
header="exact_omega omega eigenvector",
fmt=str('%.15g')
)
# Save in HDF5 format.
with h5py.File('results_py.h5', 'w') as f:
f['eigenvectors'] = eigvec
f['omega'] = omega
f['exact_omega'] = exact_omega
if __name__ == "__main__":
main()
#!/usr/bin/env python
# -*- coding: utf-8 -*-
#
# Programming Techniques for Scientific Simulations, HS 2020
# Exercise 13.2
from __future__ import (division, print_function)
import numpy as np
import scipy.stats as st
import scipy.linalg as la
import matplotlib.pyplot as plt
def get_samples(n=500):
return np.random.rand(n, 2) * 2 - 1
def get_r(samples):
return la.norm(samples, axis=1)
def calculate_pi(samples):
num_inside = np.sum(get_r(samples) < 1.)
num_total = len(samples)
return 4 * num_inside / num_total
def create_scatter_plot(samples):
fig, ax = plt.subplots()
ax.set_aspect('equal')
ax.set_xticks([-1, 0, 1])
ax.set_yticks([-1, 0, 1])
ax.set_xlabel(r'$x$')
ax.set_ylabel(r'$y$', rotation='horizontal')
ax.set_title('Points used to calculate Pi')
x, y = samples.T
ax.scatter(x, y, c=get_r(samples) < 1., cmap='bwr')
phi = np.linspace(0, 2 * np.pi, 200)
ax.plot(np.cos(phi), np.sin(phi), 'k')
fig.savefig('scatter_plot.pdf', bbox_inches='tight')
def create_contour_plot(samples):
density_func = st.gaussian_kde(samples.T, bw_method=0.2)
X, Y = np.mgrid[-1.1:1.1:200j, -1.1:1.1:200j]
positions = np.stack([X.flat, Y.flat])
Z = np.reshape(density_func(positions), X.shape)
fig, ax = plt.subplots()
ax.set_aspect('equal')
ax.set_xticks([-1, 0, 1])
ax.set_yticks([-1, 0, 1])
ax.set_xlabel(r'$x$')
ax.set_ylabel(r'$y$', rotation='horizontal')
ax.set_title('Interpolated Point Density')
heatmap = ax.contourf(X, Y, Z)
ax.contour(X, Y, Z, linewidths=0.3, colors='k')
fig.colorbar(heatmap)
fig.savefig('contour_plot.pdf', bbox_inches='tight')
def main():
np.random.seed(42)
samples = get_samples()
pi_estimate = calculate_pi(samples)
print('Pi ~= {}'.format(pi_estimate))
create_scatter_plot(samples)
create_contour_plot(samples)
if __name__ == '__main__':
main()
Markdown is supported
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