7.1 在區間[0,10]上等間距取1000個點Xi（i為下標，i=1,2,3,...,1000），並計算在這些點Xi處函式g(x)=((3x2+4x+6)sinx)/(x2+8x+6)的函式值yi(i為下標)，利用觀測點（Xi,Yi）(i=1,2,3,...,1000)，求三次樣條插值函式h(x)，並求積分g(x)從0到10和h(x)從0到10。
點選檢視程式碼
import scipy.interpolate as spi
import scipy.integrate as spi_integrate
def g(x):
return ((3 * x**2 + 4 * x + 6) * np.sin(x)) / (x**2 + 8 * x + 6)
x_values = np.linspace(0, 10, 1000)
y_values = g(x_values)
spline = spi.CubicSpline(x_values, y_values)
def h(x):
return spline(x)
integral_g, _ = spi_integrate.quad(g, 0, 10)
x_fine = np.linspace(0, 10, 10000)
y_fine = h(x_fine)
integral_h = np.trapz(y_fine, x_fine)
print(f"積分 g(x) 從 0 到 10 的結果: {integral_g}")
print(f"積分 h(x) 從 0 到 10 的結果: {integral_h}")
print("學號後四位：3028")
print("學號：2023310143028")
7.3 已知當溫度T=[700,720,740,760,780]時，過熱蒸汽體積的變化為V=[0.0977,0.1218,0.1406,0.1551,0.1664]，分別採用線性插值和三次樣條插值求解T=750,770時的體積變化，並在一個圖形介面中畫出線性插值函式和三次樣條插值函式.
點選檢視程式碼
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d, CubicSpline
T = np.array([700, 720, 740, 760, 780])
V = np.array([0.0977, 0.1218, 0.1406, 0.1551, 0.1664])
T_interp = np.array([750, 770])
f_linear = interp1d(T, V, kind='linear')
V_linear_interp = f_linear(T_interp)
cs = CubicSpline(T, V)
V_cubic_interp = cs(T_interp)
print(f"線性插值結果: T={T_interp} 對應的 V={V_linear_interp}")
print(f"三次樣條插值結果: T={T_interp} 對應的 V={V_cubic_interp}")
print("學號後四位：3004")
x = np.linspace(700, 780, 400)
plt.figure(figsize=(10, 6))
plt.plot(T, V, 'o', label='原始資料點')
plt.plot(x, f_linear(x), '-', label='線性插值')
plt.plot(x, cs(x), '--', label='三次樣條插值')
plt.scatter(T_interp, V_linear_interp, color='red', label='線性插值點')
plt.scatter(T_interp, V_cubic_interp, color='green', label='三次樣條插值點')
plt.xlabel('溫度 T')
plt.ylabel('體積 V')
plt.title('過熱蒸汽體積隨溫度變化的插值')
plt.legend()
plt.grid(True)
plt.show()
print("學號：2023310143028")
7.4 考慮矩形區域x∈[-3,3]，y∈[-4,4]上的函式f(x,y)=(x2-2x)e(-x2-y2-xy)，利用隨機生成函式uniform隨機生成該矩形內的散亂點，然後利用griddata函式進行插值處理，並作出插值結果圖形.
點選檢視程式碼
import matplotlib.pyplot as plt
from scipy.interpolate import griddata
def f(x, y):
return (x2 - 2*x) * np.exp(-x2 - y**2 - x*y)
x_min, x_max = -3, 3
y_min, y_max = -4, 4
num_points = 1000
x_random = np.random.uniform(x_min, x_max, num_points)
y_random = np.random.uniform(y_min, y_max, num_points)
z_random = f(x_random, y_random)
grid_x, grid_y = np.mgrid[x_min:x_max:100j, y_min:y_max:100j]
grid_z = griddata((x_random, y_random), z_random, (grid_x, grid_y), method='cubic')
plt.figure(figsize=(10, 8))
plt.contourf(grid_x, grid_y, grid_z, levels=50, cmap='viridis') # 使用等高線圖填充
plt.colorbar(label='f(x, y)')
plt.scatter(x_random, y_random, c='red', s=10, label='隨機散亂點') # 繪製隨機散亂點
plt.title('函式f(x, y)的插值結果')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.grid(True)
plt.show()
print("學號：2023310143028")
7.7 g(x)=(10a)/(10b+(a-10b)e^(asinx))，取a=1.1，b=0.01，計算x=1,2,...,20時，g(x)對應的函式值，把這樣得到的資料作為模擬觀測值，記作(xi,yi)（其中i為x，y的下標），i=1,2,...,20。用curve_fit擬合函式g(x)。用leastsq擬合函式g(x)。用least_squares擬合函式g(x)
點選檢視程式碼
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit, leastsq, least_squares
from scipy.constants import e # Note: 'e' is imported but not used in the code
def g(x, a, b):
return (10 * a) / (10 * b + (a - 10 * b) * np.exp(a * np.sin(x)))
a = 1.1
b = 0.01
x_values = np.arange(1, 21)
y_values = g(x_values, a, b)
for i, (xi, yi) in enumerate(zip(x_values, y_values), start=1):
print(f"({xi}, {yi:.6f})")
popt_curve_fit, pcov_curve_fit = curve_fit(g, x_values, y_values, p0=[a, b])
y_fit_curve_fit = g(x_values, *popt_curve_fit)
def func_leastsq(params, x, y):
return y - g(x, *params)
popt_leastsq = leastsq(func_leastsq, [a, b], args=(x_values, y_values))[0]
y_fit_leastsq = g(x_values, *popt_leastsq)
popt_least_squares = least_squares(func_leastsq, [a, b], args=(x_values, y_values)).x
y_fit_least_squares = g(x_values, *popt_least_squares)
print("\ncurve_fit parameters:", popt_curve_fit)
print("leastsq parameters:", popt_leastsq)
print("least_squares parameters:", popt_least_squares)
plt.figure(figsize=(10, 6))
plt.scatter(x_values, y_values, label='Simulated data', color='red')
plt.plot(x_values, y_fit_curve_fit, label='curve_fit', linestyle='-')
plt.plot(x_values, y_fit_leastsq, label='leastsq', linestyle='--')
plt.plot(x_values, y_fit_least_squares, label='least_squares', linestyle=':') # 修正線型
plt.xlabel('x')
plt.ylabel('g(x)')
plt.legend()
plt.title('Fitting of g(x) using curve_fit, leastsq, and least_squares')
plt.grid(True)
plt.show()
print("學號：2023310143028")
7.10 已知一組觀測資料，如表中7.17.excel（表中第一列為x的值，第二列為y的值）。試用插值方法繪製出x屬於-2到4.9區間內的曲線，並比較各種插值演算法的優劣。試用多項式擬合表中資料，選擇一個能較好擬合資料點的多項式的階式，給出相應多項式的係數和剩餘標準差。若表中資料滿足正態分佈函式，試用最小二乘法非線性擬合方法求出分佈引數的值，並利用所求引數值繪製擬合曲線，觀察擬合效果。
點選檢視程式碼
import pandas as pd
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d, PchipInterpolator, CubicSpline
from scipy.optimize import curve_fit
from scipy.stats import norm
file_path = '7.17.xlsx'
data = pd.read_excel(file_path, header=None)
x_data = data.iloc[:, 0].values
y_data = data.iloc[:, 1].values
x_interp = np.linspace(-2, 4.9, 400)
interp_functions = {
'Linear': interp1d(x_data, y_data, kind='linear', bounds_error=False, fill_value='extrapolate'),
'Cubic': interp1d(x_data, y_data, kind='cubic', bounds_error=False, fill_value='extrapolate'),
'PCHIP': PchipInterpolator(x_data, y_data, extrapolate=True),
'CubicSpline': CubicSpline(x_data, y_data, bc_type='natural', extrapolate=True)
}
y_interps = {name: func(x_interp) for name, func in interp_functions.items()}
plt.figure(figsize=(10, 6))
for name, y_interp in y_interps.items():
plt.plot(x_interp, y_interp, label=name)
plt.plot(x_data, y_data, 'o', label='Original Data')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.title('Interpolation Methods Comparison')
plt.grid(True)
plt.show()
degrees = range(1, 6)
coeffs = {}
residuals = {}
for degree in degrees:
coefficients, residual = np.polyfit(x_data, y_data, degree, cov=True)
residual_std = np.sqrt(residual[0, 0])
coeffs[degree] = coefficients
residuals[degree] = residual_std
for degree, coeffs_val in coeffs.items():
print(f"Degree {degree} Polynomial Coefficients: {coeffs_val}")
print(f"Residual Standard Deviation: {residuals[degree]:.4f}")
best_degree = min(residuals, key=residuals.get)
print(f"Best fitting polynomial degree: {best_degree}")
best_poly = np.poly1d(coeffs[best_degree])
y_poly_fit = best_poly(x_interp)
plt.figure(figsize=(10, 6))
plt.plot(x_interp, y_poly_fit, label=f'{best_degree}-degree Polynomial Fit')
plt.plot(x_data, y_data, 'o', label='Original Data')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.title('Polynomial Fit')
plt.grid(True)
plt.show()
def normal_dist(x, mu, sigma):
return norm.pdf(x, mu, sigma)
params, params_covariance = curve_fit(normal_dist, x_data, y_data, p0=[np.mean(x_data), np.std(x_data)])
mu, sigma = params
x_normal_fit = np.linspace(min(x_data), max(x_data), 400)
y_normal_fit = normal_dist(x_normal_fit, mu, sigma)
plt.figure(figsize=(10, 6))
plt.plot(x_normal_fit, y_normal_fit, label=f'Normal Distribution Fit\nmu={mu:.2f}, sigma={sigma:.2f}')
plt.plot(x_data, y_data, 'o', label='Original Data')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.title('Normal Distribution Fit')
plt.grid(True)
plt.show()
print("學號：2023310143028")
