# Untitled

a guest
Sep 29th, 2018
826
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. def F(C):
2.     x2 = x**2
3.     C5_2, C3_2, C1_2 = C
4.     y = (((C5_2 * x2) + C3_2) * x2 + C1_2) * x
5.     y *= 2        # scaling
6.     return y
7.
8. def Hrms(C):
9.     return 1E+08 * ((F(C) - truth)**2).sum()
10.
11. def Hmax(C):
12.     return 1E+06 * np.abs((F(C)-truth)).max()
13.
14. import numpy as np
15. from scipy.special import factorial
16. from scipy.optimize import minimize
17. import matplotlib.pyplot as plt
18.
19. quarterpi, halfpi, pi, twopi = [f*np.pi for f in (0.25, 0.5, 1, 2)]
20.
21. C5_2 =  0.03635510    # per https://space.stackexchange.com/a/30955/12102
22. C3_2 = -0.32161470
23. C1_2 =  0.78531340
24.
25. C0   = np.array([C5_2, C3_2, C1_2])
26.
27. Cnom = np.array([halfpi**5 /(2.*factorial(5)),
28.                  halfpi**3 /(2.*factorial(3)),
29.                  halfpi**1 /(2.*factorial(1))]) # if expanding about 0
30.
31. angle = np.linspace(0, halfpi, 1001)
32. truth = np.sin(angle)
33. x     = angle/halfpi   # scaling
34.
35. Cst = np.array([0.1, 0.0, 0.5])  # arbitrary starting point
36.
37. opts = {'xtol': 1E-12, 'disp': False, 'maxiter':1000, 'maxfev':2000}
38.
39. resrms = minimize(Hrms, Cst, method='Nelder-Mead', options=opts)
40. Crms   = resrms['x']
41.
42. resmax = minimize(Hmax, Cst, method='Nelder-Mead', options=opts)
43. Cmax   = resmax['x']
44.
45. np.set_printoptions(precision=6)
46.
47. print ("C0:   {}, Hmax(C0):   {:0.2f}, Hmrms(C0):   {:0.2f}".format(C0, Hmax(C0), Hrms(C0)))
48. print ("Cnom: {}, Hmax(Cnom): {:0.2f}, Hmrms(Cnom): {:0.2f}".format(Cnom, Hmax(Cnom), Hrms(Cnom)))
49. print ("Crms: {}, Hmax(Crms): {:0.2f}, Hmrms(Crms): {:0.2f}".format(Crms, Hmax(Crms), Hrms(Crms)))
50. print ("Cmax: {}, Hmax(Cmax): {:0.2f}, Hmrms(Cmax): {:0.2f}".format(Cmax, Hmax(Cmax), Hrms(Cmax)))