import os
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
os.chdir(r'D:\data\scripts\wordpress\ex55')
os.getcwd()
x0,A,gamma=12,3,5
n=200
x=np.linspace(1,20,n)
yexact=A*gamma**2/(gamma**2+(x-x0)**2)
#Learning Scientific Programming with Python
# Add some noise with a sigma of 0.5 apart from a particularly noisy region
# near x0 where sigma is 3
sigma=np.ones(n)*0.5
sigma[np.abs(x-x0+1)<1]=3
noise=np.random.randn(n)*sigma
y=yexact+noise

def f(x,x0,A,gamma):
    #The Lorentzian entered at x0 with amplitude A and HWHM gamma
    return A*gamma**2/(gamma**2+(x-x0)**2)

def rms(y,yfit):
    return np.sqrt(np.sum((y-yfit)**2))

# Unweighted fit
p0=10,4,2
popt,pcov=curve_fit(f,x,y,p0)
yfit=f(x,*popt)
print('Unweighted fit parameters:',popt)
print('Covarian cematrix:')
print(pcov)
print('rms error in fit:',rms(yexact,yfit))
print()
# Weighted fit
popt2,pcov2=curve_fit(f,x,y,p0,sigma=sigma,absolute_sigma=True)
yfit2=f(x,*popt2)
print('Weighted fit parameters:',popt2)
print('Covariancematrix:')
print(pcov2)
print('rms error in fit:',rms(yexact,yfit2))
plt.plot(x,yexact,label='Exact')
plt.plot(x,y,'o',label='Noisy data')
plt.plot(x,yfit,label='Unweighted fit')
plt.plot(x,yfit2,label='Weighted fit')
plt.ylim(-1,4)
plt.legend(loc='lower center')

plt.tight_layout()
plt.savefig("example55.png", dpi=100)
plt.show()
plt.close()



example55

Discover more from Tips and Hints for Aerospace Engineers

Subscribe now to keep reading and get access to the full archive.

Continue reading