第七章習題

學號後四位:3018
7.3:

點選檢視程式碼

import numpy as np
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])

# 線性插值
linear_interp = interp1d(T, V)
V_linear_750 = linear_interp(750)
V_linear_770 = linear_interp(770)

# 三次樣條插值
cubic_interp = CubicSpline(T, V)
V_cubic_750 = cubic_interp(750)
V_cubic_770 = cubic_interp(770)

print(f"線性插值：T = 750 時，V = {V_linear_750}；T = 770 時，V = {V_linear_770}")
print(f"三次樣條插值：T = 750 時，V = {V_cubic_750}；T = 770 時，V = {V_cubic_770}")

# 繪製圖形
T_new = np.linspace(700, 780, 100)
V_linear_new = linear_interp(T_new)
V_cubic_new = cubic_interp(T_new)

plt.plot(T, V, 'o', label='原始資料')
plt.plot(T_new, V_linear_new, label='線性插值')
plt.plot(T_new, V_cubic_new, label='三次樣條插值')
plt.xlabel('溫度 T')
plt.ylabel('體積變化 V')
plt.legend()
plt.show()
print("xuehao3018")

7.4:

點選檢視程式碼

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import griddata

# 定義函式
def f(x, y):
    return (x**2 - 2*x) * np.exp(-y**2/2)

# 生成隨機散亂點
np.random.seed(0)
n_points = 100
x_random = np.random.uniform(-3, 3, n_points)
y_random = np.random.uniform(-4, 4, n_points)
z_random = f(x_random, y_random)

# 定義插值網格
grid_x, grid_y = np.mgrid[-3:3:100j, -4:4:100j]

# 進行插值
z_interpolated = griddata((x_random, y_random), z_random, (grid_x, grid_y), method='cubic')

# 繪製結果
plt.figure(figsize=(10, 5))
plt.subplot(121)
plt.scatter(x_random, y_random, c=z_random, cmap='viridis')
plt.title('Scattered Points')
plt.xlabel('x')
plt.ylabel('y')

plt.subplot(122)
plt.imshow(z_interpolated.T, extent=[-3, 3, -4, 4], origin='lower', cmap='viridis')
plt.title('Interpolated Result')
plt.xlabel('x')
plt.ylabel('y')

plt.show()
print("xuehao3018")

7.7：

點選檢視程式碼

import numpy as np
from scipy.optimize import curve_fit, least_squares
import matplotlib.pyplot as plt

# 定義函式 g(x)
def g(x, j0, q, r):
    return j0 * (1.1 ** x) + (q - 1.1 * j0) * r * x

# 生成模擬觀測值
a = 1.1
b = 0.01
x_data = np.arange(1, 21)
y_data = [g(x, a, a + b, b) for x in x_data]

# 1. 用 curve_fit 擬合函式 g(x)
popt_curve_fit, _ = curve_fit(g, x_data, y_data)
y_curve_fit = g(x_data, *popt_curve_fit)

# 2. 用 leastsq 擬合函式 g(x)（注意：leastsq 在較新版本的 scipy 中已被棄用，這裡只是為了演示）
from scipy.optimize import leastsq

def residuals(params, x, y):
    return y - g(x, *params)

initial_guess = [1, 1, 1]
popt_leastsq, _ = leastsq(residuals, initial_guess, args=(x_data, y_data))
y_leastsq = g(x_data, *popt_leastsq)

# 3. 用 least_squares 擬合函式 g(x)
def fun(params):
    return g(x_data, *params) - y_data

initial_guess_ls = [1, 1, 1]
res_ls = least_squares(fun, initial_guess_ls)
popt_least_squares = res_ls.x
y_least_squares = g(x_data, *popt_least_squares)

# 繪製結果
plt.plot(x_data, y_data, 'o', label='模擬觀測值')
plt.plot(x_data, y_curve_fit, label='curve_fit')
plt.plot(x_data, y_leastsq, label='leastsq')
plt.plot(x_data, y_least_squares, label='least_squares')
plt.legend()
plt.xlabel('x')
plt.ylabel('g(x)')
plt.show()
print("xuehao3018")

7.10:

點選檢視程式碼

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d
from scipy.optimize import curve_fit

# 資料
x_data = np.array([-2, -1.7, -1.4, -1.1, -0.8, -0.5, -0.2, 0.1, 0.4, 0.7, 1, 1.3, 1.6, 1.9, 2.2, 2.5, 2.8, 3.1, 3.4, 3.7, 4, 4.3, 4.6, 4.9])
y_data = np.array([0.1029, 0.1174, 0.1316, 0.1448, 0.1566, 0.1662, 0.1733, 0.1775, 0.1785, 0.1764, 0.1711, 0.1630, 0.1526, 0.1402, 0.1266, 0.1122, 0.0977, 0.0835, 0.0702, 0.0588, 0.0479, 0.0373, 0.0291, 0.0224])

# 1. 插值方法繪製曲線並比較
# 線性插值
linear_interp = interp1d(x_data, y_data, kind='linear')
x_linear = np.linspace(-2, 4.9, 100)
y_linear = linear_interp(x_linear)

# 三次樣條插值
cubic_interp = interp1d(x_data, y_data, kind='cubic')
x_cubic = np.linspace(-2, 4.9, 100)
y_cubic = cubic_interp(x_cubic)

# 繪製插值結果
plt.figure(figsize=(10, 6))
plt.plot(x_data, y_data, 'o', label='原始資料')
plt.plot(x_linear, y_linear, label='線性插值')
plt.plot(x_cubic, y_cubic, label='三次樣條插值')
plt.xlabel('x')
plt.ylabel('y')
plt.title('插值方法比較')
plt.legend()

# 比較優劣（簡單描述）
print("線性插值簡單快速，但在曲線變化劇烈處可能不夠準確。三次樣條插值通常更光滑，但計算量較大。")

# 2. 多項式擬合
def polynomial_func(x, *coeffs):
    y = 0
    for i, c in enumerate(coeffs):
        y += c * x**i
    return y

degrees = [2, 3, 4, 5]
best_degree = None
best_residuals = None
best_coeffs = None

for degree in degrees:
    initial_coeffs = np.ones(degree + 1)
    popt, _ = curve_fit(polynomial_func, x_data, y_data, p0=initial_coeffs)
    residuals = y_data - polynomial_func(x_data, *popt)
    if best_residuals is None or np.sum(residuals**2) < best_residuals:
        best_degree = degree
        best_residuals = np.sum(residuals**2)
        best_coeffs = popt

print(f"最佳多項式階次為 {best_degree}，係數為 {best_coeffs}，剩餘標準差為 {np.sqrt(best_residuals / len(x_data))}")

# 3. 正態分佈擬合
def normal_dist_func(x, mu, sigma):
    return 1 / (sigma * np.sqrt(2 * np.pi)) * np.exp(-(x - mu)**2 / (2 * sigma**2))

popt_normal, _ = curve_fit(normal_dist_func, x_data, y_data)
mu_fit, sigma_fit = popt_normal

x_normal = np.linspace(-2, 4.9, 100)
y_normal = normal_dist_func(x_normal, mu_fit, sigma_fit)

plt.figure(figsize=(10, 6))
plt.plot(x_data, y_data, 'o', label='原始資料')
plt.plot(x_normal, y_normal, label='正態分佈擬合')
plt.xlabel('x')
plt.ylabel('y')
plt.title('正態分佈擬合')
plt.legend()

plt.show()
print("xuehao3018")