import numpy as np

# 힐버트 행렬 (악명 높은 나쁜 조건 행렬)
n = 6
H = np.array([[1/(i+j+1) for j in range(n)] for i in range(n)])
kappa = np.linalg.cond(H)
print(f"힐버트 행렬 ({n}x{n}) 조건수: {kappa:.2e}")

# 간단한 행렬
A = np.array([[2, 1], [1, 3]], dtype=float)
print(f"간단한 행렬 조건수: {np.linalg.cond(A):.2f}")

# 오차 전파 시연
b_true = np.array([3, 4], dtype=float)
x_true = np.linalg.solve(A, b_true)

# 약간의 오차 추가
b_perturbed = b_true + np.array([0.01, 0.0])
x_perturbed = np.linalg.solve(A, b_perturbed)

print(f"\n입력 오차: {np.linalg.norm(b_perturbed - b_true):.4f}")
print(f"출력 오차: {np.linalg.norm(x_perturbed - x_true):.4f}")

def bisection(f, a, b, tol=1e-10, max_iter=100):
    """이분법으로 f(x) = 0의 근 찾기"""
    if f(a) * f(b) > 0:
        raise ValueError("f(a)와 f(b)는 반대 부호여야 합니다")

    history = []
    for i in range(max_iter):
        c = (a + b) / 2
        history.append({'iter': i+1, 'a': a, 'b': b, 'c': c,
                        'f(c)': f(c), 'error': abs(b-a)/2})

        if abs(f(c)) < tol or (b - a) / 2 < tol:
            print(f"수렴: {i+1}번 반복")
            return c, history

        if f(a) * f(c) < 0:
            b = c
        else:
            a = c

    return (a + b) / 2, history

# 예제: x^3 - 2x - 5 = 0
def equation(x):
    return x**3 - 2*x - 5

root, history = bisection(equation, 1, 3)
print(f"근: {root:.10f}")
print(f"검증: f({root:.6f}) = {equation(root):.2e}")

# 수렴 이력 출력
print("\n반복 | 근사값 | 오차")
for h in history[:8]:
    print(f"  {h['iter']:2d} | {h['c']:.8f} | {h['error']:.2e}")

def newton_raphson(f, df, x0, tol=1e-10, max_iter=50):
    """뉴턴-랩슨법"""
    x = x0
    history = [{'iter': 0, 'x': x, 'f(x)': f(x), 'error': float('inf')}]

    for i in range(max_iter):
        fx = f(x)
        dfx = df(x)

        if abs(dfx) < 1e-15:
            print("경고: 도함수가 0에 가까움, 발산 위험")
            break

        x_new = x - fx / dfx
        error = abs(x_new - x)
        history.append({'iter': i+1, 'x': x_new,
                        'f(x)': f(x_new), 'error': error})

        if error < tol:
            print(f"수렴: {i+1}번 반복")
            return x_new, history

        x = x_new

    return x, history

# f(x) = x^3 - 2x - 5, f'(x) = 3x^2 - 2
def f(x):
    return x**3 - 2*x - 5

def df(x):
    return 3*x**2 - 2

root, history = newton_raphson(f, df, x0=2.0)
print(f"뉴턴-랩슨 근: {root:.10f}")

print("\n반복 | 근사값 | 오차")
for h in history:
    err_str = f"{h['error']:.2e}" if h['error'] != float('inf') else "   -"
    print(f"  {h['iter']:2d} | {h['x']:.10f} | {err_str}")

import numpy as np
import matplotlib.pyplot as plt

# 이분법 vs 뉴턴-랩슨 수렴 속도 비교
true_root = 2.0945514815423265  # x^3 - 2x - 5 = 0의 실제 근

_, hist_bisect = bisection(f, 1, 3, tol=1e-12, max_iter=50)
_, hist_newton = newton_raphson(f, df, x0=2.0, tol=1e-12, max_iter=20)

errors_b = [abs(h['c'] - true_root) for h in hist_bisect]
errors_n = [abs(h['x'] - true_root) for h in hist_newton if h['error'] != float('inf')]

plt.figure(figsize=(10, 5))
plt.semilogy(range(1, len(errors_b)+1), errors_b, 'b-o', label='이분법 (1차 수렴)')
plt.semilogy(range(1, len(errors_n)+1), errors_n, 'r-s', label='뉴턴-랩슨 (2차 수렴)')
plt.xlabel('반복 횟수')
plt.ylabel('절대 오차 (로그 스케일)')
plt.title('이분법 vs 뉴턴-랩슨 수렴 속도')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('convergence_comparison.png', dpi=150)
plt.show()

def secant_method(f, x0, x1, tol=1e-10, max_iter=50):
    """할선법"""
    for i in range(max_iter):
        f0, f1 = f(x0), f(x1)
        if abs(f1 - f0) < 1e-15:
            break
        x2 = x1 - f1 * (x1 - x0) / (f1 - f0)
        if abs(x2 - x1) < tol:
            print(f"수렴: {i+1}번 반복")
            return x2
        x0, x1 = x1, x2
    return x1

root_s = secant_method(f, 1.0, 3.0)
print(f"할선법 근: {root_s:.10f}")

from scipy.optimize import fsolve, brentq, minimize_scalar
import numpy as np

# 단일 방정식 풀기
def equation(x):
    return x**3 - 2*x - 5

# 브렌트법 (이분법 + 할선법 + 역 이차 보간 혼합)
root = brentq(equation, 1, 3, xtol=1e-12)
print(f"brentq 근: {root:.12f}")

# fsolve (뉴턴-랩슨 기반)
root_fs = fsolve(equation, 2.0)[0]
print(f"fsolve 근: {root_fs:.12f}")

# 연립 비선형 방정식
from scipy.optimize import fsolve

def system(vars):
    x, y = vars
    eq1 = x**2 + y**2 - 4   # 원
    eq2 = x*y - 1             # 쌍곡선
    return [eq1, eq2]

solutions = []
for x0, y0 in [(1, 1), (-1, 1), (1, -1), (-1, -1)]:
    sol = fsolve(system, [x0, y0], full_output=True)
    if sol[2] == 1:  # 수렴 성공
        solutions.append(sol[0])

print("\n연립 방정식 해:")
for sol in solutions:
    print(f"  (x, y) = ({sol[0]:.6f}, {sol[1]:.6f})")

# 회로 방정식: 다이오드 특성 (비선형)
def diode_equation(V, V_s=5, R=1000, I_s=1e-12, n=1, V_T=0.02585):
    """V_s = V + R * I_s * (exp(V/(n*V_T)) - 1)"""
    I = I_s * (np.exp(V / (n * V_T)) - 1)
    return V_s - V - R * I

V_diode = brentq(diode_equation, 0, 1)
I_diode = 1e-12 * (np.exp(V_diode / 0.02585) - 1)
print(f"\n다이오드 동작점: V = {V_diode:.4f} V, I = {I_diode*1000:.4f} mA")

import numpy as np
import matplotlib.pyplot as plt

def numerical_derivative(f, x, h=1e-5, method='central'):
    """수치 미분"""
    if method == 'forward':
        return (f(x + h) - f(x)) / h
    elif method == 'backward':
        return (f(x) - f(x - h)) / h
    elif method == 'central':
        return (f(x + h) - f(x - h)) / (2 * h)

def f(x):
    return np.sin(x)

def df_true(x):
    return np.cos(x)

x0 = np.pi / 4
true_deriv = df_true(x0)

# h 크기에 따른 오차 분석
h_values = np.logspace(-16, 0, 100)
errors_fwd = []
errors_cen = []

for h in h_values:
    e_fwd = abs(numerical_derivative(f, x0, h, 'forward') - true_deriv)
    e_cen = abs(numerical_derivative(f, x0, h, 'central') - true_deriv)
    errors_fwd.append(e_fwd if e_fwd > 0 else 1e-17)
    errors_cen.append(e_cen if e_cen > 0 else 1e-17)

plt.figure(figsize=(10, 6))
plt.loglog(h_values, errors_fwd, 'b-', linewidth=2, label='전진 차분 O(h)')
plt.loglog(h_values, errors_cen, 'r-', linewidth=2, label='중앙 차분 O(h^2)')
plt.loglog(h_values, h_values, 'b--', alpha=0.5, label='O(h) 기준선')
plt.loglog(h_values, h_values**2, 'r--', alpha=0.5, label='O(h^2) 기준선')
plt.xlabel('h (스텝 크기)')
plt.ylabel('절대 오차')
plt.title('수치 미분 오차 분석 (h에 따른 오차)')
plt.legend()
plt.grid(True, which='both', alpha=0.3)
plt.tight_layout()
plt.savefig('numerical_derivative_error.png', dpi=150)
plt.show()

import numpy as np
from scipy import integrate
import matplotlib.pyplot as plt

def trapezoidal(f, a, b, n):
    """사다리꼴 법칙"""
    x = np.linspace(a, b, n+1)
    h = (b - a) / n
    y = f(x)
    return h * (y[0]/2 + np.sum(y[1:-1]) + y[-1]/2)

def simpsons(f, a, b, n):
    """심프슨 1/3 법칙 (n은 짝수)"""
    if n % 2 != 0:
        n += 1
    x = np.linspace(a, b, n+1)
    h = (b - a) / n
    y = f(x)
    return h/3 * (y[0] + 4*np.sum(y[1::2]) + 2*np.sum(y[2:-1:2]) + y[-1])

# 예제: integral of sin(x) from 0 to pi = 2
f_test = np.sin
a, b = 0, np.pi
true_val = 2.0

print("수치 적분 정확도 비교")
print(f"정확한 값: {true_val}")
print("-" * 60)
print(f"{'방법':<20} {'n':>6} {'근사값':>15} {'상대 오차':>12}")
print("-" * 60)

for n in [4, 8, 16, 32]:
    trap = trapezoidal(f_test, a, b, n)
    simp = simpsons(f_test, a, b, n)
    print(f"{'사다리꼴':<20} {n:>6} {trap:>15.10f} {abs(trap-true_val)/true_val:>12.2e}")
    print(f"{'심프슨':<20} {n:>6} {simp:>15.10f} {abs(simp-true_val)/true_val:>12.2e}")

# scipy.integrate.quad - 적응형 구적법
result, error = integrate.quad(f_test, a, b)
print(f"{'scipy.quad':<20} {'적응형':>6} {result:>15.10f} {'<'+str(error):>12}")

# 어려운 적분: int_0^inf exp(-x^2) dx = sqrt(pi)/2
def gaussian(x):
    return np.exp(-x**2)

result_gauss, _ = integrate.quad(gaussian, 0, np.inf)
print(f"\n가우스 적분: {result_gauss:.10f}")
print(f"정확한 값 (sqrt(pi)/2): {np.sqrt(np.pi)/2:.10f}")

# 이중 적분: int_0^1 int_0^1 x*y dx dy = 0.25
def f2d(y, x):
    return x * y

result_2d, _ = integrate.dblquad(f2d, 0, 1, 0, 1)
print(f"\n이중 적분: {result_2d:.10f} (정확: 0.25)")

import numpy as np
import matplotlib.pyplot as plt

def euler_method(f, x0, y0, h, x_end):
    """오일러법"""
    x_vals = [x0]
    y_vals = [y0]
    x = x0
    y = y0

    while x < x_end - 1e-10:
        h_actual = min(h, x_end - x)
        y = y + h_actual * f(x, y)
        x = x + h_actual
        x_vals.append(x)
        y_vals.append(y)

    return np.array(x_vals), np.array(y_vals)

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

def rk4_method(f, x0, y0, h, x_end):
    """4차 룽게-쿠타법"""
    x_vals = [x0]
    y_vals = [y0]
    x = x0
    y = y0

    while x < x_end - 1e-10:
        h_actual = min(h, x_end - x)
        k1 = h_actual * f(x, y)
        k2 = h_actual * f(x + h_actual/2, y + k1/2)
        k3 = h_actual * f(x + h_actual/2, y + k2/2)
        k4 = h_actual * f(x + h_actual, y + k3)
        y = y + (k1 + 2*k2 + 2*k3 + k4) / 6
        x = x + h_actual
        x_vals.append(x)
        y_vals.append(y)

    return np.array(x_vals), np.array(y_vals)

# 완전한 RLC 회로 시뮬레이션 - 여러 방법 비교
R, L, C = 50, 0.2, 100e-6
Vs = 100

def rlc_ode(t, state):
    q, dq = state
    d2q = (Vs - R * dq - q/C) / L
    return [dq, d2q]

def rlc_scalar(t, y):
    """연립 ODE를 단일 함수로"""
    return rlc_ode(t, y)

# 초기조건
y0 = [0.0, 0.0]
t_end = 0.02
h_values = [0.001, 0.0005, 0.0001]

# scipy solve_ivp (기준 해)
t_ref = np.linspace(0, t_end, 5000)
sol_ref = solve_ivp(rlc_ode, [0, t_end], y0,
                    t_eval=t_ref, method='RK45',
                    rtol=1e-10, atol=1e-12)

v_ref = sol_ref.y[0] / C  # 커패시터 전압

# RK4로 계산
def rlc_2nd_order(t, y):
    """연립 1차 시스템"""
    q, dq = y[0], y[1]
    return [dq, (Vs - R*dq - q/C) / L]

fig, axes = plt.subplots(2, 1, figsize=(12, 10))

axes[0].plot(sol_ref.t * 1000, v_ref, 'k-', linewidth=2.5, label='기준 (RK45 고정밀)')

for h in h_values:
    n_steps = int(t_end / h)
    t_rk4 = np.linspace(0, t_end, n_steps + 1)

    # 수동 RK4
    y = np.array(y0, dtype=float)
    y_hist = [y.copy()]
    for i in range(n_steps):
        t_i = t_rk4[i]
        k1 = h * np.array(rlc_2nd_order(t_i, y))
        k2 = h * np.array(rlc_2nd_order(t_i + h/2, y + k1/2))
        k3 = h * np.array(rlc_2nd_order(t_i + h/2, y + k2/2))
        k4 = h * np.array(rlc_2nd_order(t_i + h, y + k3))
        y = y + (k1 + 2*k2 + 2*k3 + k4) / 6
        y_hist.append(y.copy())

    y_arr = np.array(y_hist)
    v_rk4 = y_arr[:, 0] / C
    axes[0].plot(t_rk4 * 1000, v_rk4, '--',
                 label=f'RK4 (h={h*1000:.1f}ms)', alpha=0.8)

axes[0].set_xlabel('시간 (ms)')
axes[0].set_ylabel('커패시터 전압 (V)')
axes[0].set_title('RLC 직렬 회로 - RK4 수치 풀이')
axes[0].legend()
axes[0].grid(True, alpha=0.3)
axes[0].axhline(y=Vs, color='gray', linestyle=':', alpha=0.7, label='정상 상태')

# 전류 플롯
i_ref = sol_ref.y[1]
axes[1].plot(sol_ref.t * 1000, i_ref * 1000, 'k-', linewidth=2, label='전류 i(t)')
axes[1].set_xlabel('시간 (ms)')
axes[1].set_ylabel('전류 (mA)')
axes[1].set_title('RLC 회로 전류')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('rlc_numerical.png', dpi=150)
plt.show()

# 감쇠 특성 확인
omega_0 = 1 / np.sqrt(L * C)
alpha = R / (2 * L)
omega_d = np.sqrt(max(omega_0**2 - alpha**2, 0))
print(f"공진 주파수: {omega_0/(2*np.pi):.1f} Hz")
print(f"감쇠 계수: {alpha:.1f} s^-1")
if omega_0 > alpha:
    print(f"부족감쇠, 감쇠 진동 주파수: {omega_d/(2*np.pi):.1f} Hz")
elif omega_0 == alpha:
    print("임계감쇠")
else:
    print("과감쇠")

from scipy.integrate import solve_ivp
import numpy as np

# 강성 방정식 예시: Robertson 화학 반응
def robertson(t, y):
    """
    dy0/dt = -0.04*y0 + 1e4*y1*y2
    dy1/dt = 0.04*y0 - 1e4*y1*y2 - 3e7*y1^2
    dy2/dt = 3e7*y1^2
    """
    dy0 = -0.04*y[0] + 1e4*y[1]*y[2]
    dy1 = 0.04*y[0] - 1e4*y[1]*y[2] - 3e7*y[1]**2
    dy2 = 3e7*y[1]**2
    return [dy0, dy1, dy2]

y0 = [1, 0, 0]
t_span = (0, 1e11)
t_eval = np.logspace(-6, 11, 1000)

# 비강성 방법 (RK45): 너무 느림
# 강성 방법 (Radau, BDF): 효율적

sol = solve_ivp(robertson, t_span, y0,
                method='Radau',  # 강성 방법
                t_eval=t_eval,
                rtol=1e-6, atol=1e-9)

print(f"함수 평가 횟수: {sol.nfev}")
print(f"야코비안 평가: {sol.njev}")

import numpy as np

def gaussian_elimination_pivoting(A, b):
    """부분 피벗팅을 포함한 가우스 소거법"""
    n = len(b)
    A = np.array(A, dtype=float)
    b = np.array(b, dtype=float)

    # 전진 소거
    for i in range(n):
        # 부분 피벗팅: i열에서 최대 원소 찾기
        max_row = i + np.argmax(abs(A[i:, i]))
        A[[i, max_row]] = A[[max_row, i]]
        b[[i, max_row]] = b[[max_row, i]]

        if abs(A[i, i]) < 1e-15:
            raise ValueError("특이 행렬!")

        for j in range(i+1, n):
            factor = A[j, i] / A[i, i]
            A[j, i:] -= factor * A[i, i:]
            b[j] -= factor * b[i]

    # 후진 대입
    x = np.zeros(n)
    for i in range(n-1, -1, -1):
        x[i] = (b[i] - np.dot(A[i, i+1:], x[i+1:])) / A[i, i]

    return x

# 예제: 3개 루프 회로망 (KVL)
A = np.array([
    [10, -3, -1],
    [-3, 12, -2],
    [-1, -2,  8]
], dtype=float)
b = np.array([100, 50, -30], dtype=float)

x_gauss = gaussian_elimination_pivoting(A, b)
x_numpy = np.linalg.solve(A, b)

print("회로 루프 전류 (가우스 소거법):")
for i, v in enumerate(x_gauss):
    print(f"  I{i+1} = {v:.6f} A")

print(f"\nnumpy 결과와 오차: {np.max(abs(x_gauss - x_numpy)):.2e}")

from scipy.linalg import lu, lu_factor, lu_solve
import numpy as np

A = np.array([[2, 1, -1], [-3, -1, 2], [-2, 1, 2]], dtype=float)
b = np.array([8, -11, -3], dtype=float)

# LU 분해
P, L, U = lu(A)
print("L 행렬:")
print(np.round(L, 4))
print("U 행렬:")
print(np.round(U, 4))

# LU 분해로 풀기
lu_and_piv = lu_factor(A)
x = lu_solve(lu_and_piv, b)
print(f"\n해: {x}")

# 여러 우변 벡터에 대해 효율적으로 풀기
for b_new in [b, b+1, b*2]:
    x_new = lu_solve(lu_and_piv, b_new)  # L, U 재사용

import numpy as np

def gauss_seidel(A, b, x0=None, tol=1e-10, max_iter=1000):
    """가우스-자이델 반복법"""
    n = len(b)
    x = np.zeros(n) if x0 is None else x0.copy()

    for iteration in range(max_iter):
        x_old = x.copy()
        for i in range(n):
            sum1 = np.dot(A[i, :i], x[:i])
            sum2 = np.dot(A[i, i+1:], x_old[i+1:])
            x[i] = (b[i] - sum1 - sum2) / A[i, i]

        residual = np.linalg.norm(x - x_old)
        if residual < tol:
            print(f"수렴: {iteration+1}번 반복, 잔차: {residual:.2e}")
            return x

    print(f"최대 반복 {max_iter}에서 종료, 잔차: {residual:.2e}")
    return x

# 대각 우세 행렬 예제 (노드 전압 해석)
A = np.array([
    [ 4, -1, -1,  0],
    [-1,  4,  0, -1],
    [-1,  0,  4, -1],
    [ 0, -1, -1,  4]
], dtype=float)
b = np.array([5, 10, 15, 20], dtype=float)

x_gs = gauss_seidel(A, b)
x_exact = np.linalg.solve(A, b)
print(f"가우스-자이델: {x_gs}")
print(f"정확한 해: {x_exact}")
print(f"최대 오차: {np.max(abs(x_gs - x_exact)):.2e}")

import numpy as np

A = np.array([[3, -0.1, -0.2],
              [0.1, 7, -0.3],
              [0.3, -0.2, 10]], dtype=float)
b = np.array([7.85, -19.3, 71.4], dtype=float)

# 직접 풀기
x = np.linalg.solve(A, b)
print(f"해: {x}")

# 고유값
eigenvalues, eigenvectors = np.linalg.eig(A)
print(f"\n고유값: {eigenvalues}")

# 역행렬 (비추천 - 수치 오차 큼)
A_inv = np.linalg.inv(A)
print(f"역행렬 검증 (A*A_inv - I 최대값): {np.max(abs(A @ A_inv - np.eye(3))):.2e}")

# 조건수와 랭크
print(f"조건수: {np.linalg.cond(A):.2f}")
print(f"랭크: {np.linalg.matrix_rank(A)}")

# 최소자승 풀이 (과결정 시스템)
A_over = np.vstack([A, [1, 1, 1]])
b_over = np.append(b, 10)
x_ls, residuals, rank, sv = np.linalg.lstsq(A_over, b_over, rcond=None)
print(f"\n최소자승 해: {x_ls}")

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import lagrange, CubicSpline, interp1d

# 측정 데이터 (온도 vs 저항)
x_data = np.array([0, 20, 40, 60, 80, 100], dtype=float)  # 섭씨
y_data = np.array([100, 108, 116, 125, 134, 143], dtype=float)  # Ohm

x_fine = np.linspace(0, 100, 300)

# 라그랑주 보간
poly = lagrange(x_data, y_data)
y_lagrange = poly(x_fine)

# 3차 스플라인
cs = CubicSpline(x_data, y_data)
y_spline = cs(x_fine)

# 선형 보간
f_linear = interp1d(x_data, y_data)
y_linear = f_linear(x_fine)

plt.figure(figsize=(10, 6))
plt.scatter(x_data, y_data, s=100, zorder=5, label='측정 데이터', color='black')
plt.plot(x_fine, y_linear, 'g--', linewidth=1.5, label='선형 보간')
plt.plot(x_fine, y_lagrange, 'r-', linewidth=2, label='라그랑주 보간')
plt.plot(x_fine, y_spline, 'b-', linewidth=2, label='3차 스플라인')
plt.xlabel('온도 (C)')
plt.ylabel('저항 (Ohm)')
plt.title('온도-저항 보간 비교')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('interpolation.png', dpi=150)
plt.show()

# 50도에서의 저항 예측
print(f"50도에서 저항 예측:")
print(f"  선형: {f_linear(50):.2f} Ohm")
print(f"  스플라인: {cs(50):.2f} Ohm")

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

# 노이즈가 있는 지수 감쇠 데이터
np.random.seed(42)
t_data = np.linspace(0, 5, 30)
y_true = 3 * np.exp(-0.7 * t_data)
y_noisy = y_true + 0.2 * np.random.randn(len(t_data))

# 다항식 최소자승 피팅
degrees = [1, 2, 3, 5]
t_fine = np.linspace(0, 5, 200)

plt.figure(figsize=(12, 8))
plt.scatter(t_data, y_noisy, s=50, color='black', zorder=5, label='측정 데이터')
plt.plot(t_fine, 3*np.exp(-0.7*t_fine), 'k-', linewidth=2, label='참 신호')

for deg in degrees:
    coeffs = np.polyfit(t_data, y_noisy, deg)
    y_fit = np.polyval(coeffs, t_fine)
    plt.plot(t_fine, y_fit, '--', linewidth=1.5, label=f'{deg}차 다항식')

plt.ylim([-1, 4])
plt.xlabel('시간 (s)')
plt.ylabel('진폭')
plt.title('최소자승 다항식 피팅')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('polynomial_fit.png', dpi=150)
plt.show()

# 비선형 최소자승: 지수 함수 피팅
def exp_model(t, A, k):
    return A * np.exp(-k * t)

popt, pcov = curve_fit(exp_model, t_data, y_noisy, p0=[3, 0.5])
A_fit, k_fit = popt
A_std, k_std = np.sqrt(np.diag(pcov))

print(f"지수 피팅 결과:")
print(f"  A = {A_fit:.4f} +/- {A_std:.4f} (참값: 3)")
print(f"  k = {k_fit:.4f} +/- {k_std:.4f} (참값: 0.7)")

# 잔차 및 R^2
y_pred = exp_model(t_data, *popt)
ss_res = np.sum((y_noisy - y_pred)**2)
ss_tot = np.sum((y_noisy - np.mean(y_noisy))**2)
r_squared = 1 - ss_res / ss_tot
print(f"  R^2 = {r_squared:.6f}")

import numpy as np

def power_method(A, tol=1e-10, max_iter=1000):
    """거듭제곱법 (최대 고유값)"""
    n = A.shape[0]
    q = np.ones(n) / np.sqrt(n)  # 초기 벡터

    for k in range(max_iter):
        z = A @ q
        eigenvalue = np.dot(q, z)  # 레일리 몫
        q_new = z / np.linalg.norm(z)

        if np.linalg.norm(q_new - q) < tol or np.linalg.norm(q_new + q) < tol:
            print(f"수렴: {k+1}번 반복")
            return eigenvalue, q_new

        q = q_new

    return eigenvalue, q

# 진동 시스템의 고유 주파수 분석
# 이중 진자 질량-스프링 시스템의 강성 행렬
k = 1000  # N/m
m = 1     # kg

K = np.array([[2*k, -k], [-k, k]], dtype=float)
M = np.array([[m, 0], [0, m]], dtype=float)

# 일반화 고유값 문제: K*v = lambda*M*v
# 변환: M^(-1/2) * K * M^(-1/2) * u = lambda * u
M_inv_sqrt = np.diag(1/np.sqrt(np.diag(M)))
A_sym = M_inv_sqrt @ K @ M_inv_sqrt

lam_max, v_max = power_method(A_sym)
omega_max = np.sqrt(lam_max)

# 모든 고유값 (numpy)
eigenvalues_all = np.linalg.eigvalsh(A_sym)
omega_all = np.sqrt(eigenvalues_all)

print("진동 시스템 고유 주파수:")
for i, omega in enumerate(omega_all):
    print(f"  모드 {i+1}: {omega/(2*np.pi):.3f} Hz")

print(f"\n거듭제곱법 최대 고유 주파수: {omega_max/(2*np.pi):.3f} Hz")

import numpy as np

def qr_algorithm(A, max_iter=1000, tol=1e-10):
    """기본 QR 알고리즘"""
    A_k = A.copy().astype(float)
    n = A.shape[0]

    for k in range(max_iter):
        Q, R = np.linalg.qr(A_k)
        A_k_new = R @ Q

        # 하삼각 원소들의 크기 확인
        off_diag_norm = np.max(np.abs(np.tril(A_k_new, -1)))
        if off_diag_norm < tol:
            print(f"QR 수렴: {k+1}번 반복")
            break

        A_k = A_k_new

    return np.diag(A_k)

# 대칭 행렬 (실수 고유값)
A = np.array([
    [4, 1, 0],
    [1, 3, 1],
    [0, 1, 2]
], dtype=float)

eig_qr = qr_algorithm(A)
eig_numpy = np.linalg.eigvalsh(A)

print("QR 알고리즘 고유값:", np.sort(eig_qr))
print("numpy 고유값:    ", np.sort(eig_numpy))

import numpy as np
import matplotlib.pyplot as plt
from scipy import linalg

def vibration_analysis(K, M, n_modes=None):
    """
    일반화 고유값 문제: K*v = omega^2 * M * v
    K: 강성 행렬, M: 질량 행렬
    """
    # scipy linalg.eigh로 일반화 고유값 풀기
    eigenvalues, eigenvectors = linalg.eigh(K, M)

    omega = np.sqrt(np.abs(eigenvalues))  # 각 주파수 (rad/s)
    freq = omega / (2 * np.pi)           # 주파수 (Hz)

    if n_modes is None:
        n_modes = len(eigenvalues)

    print("고유 주파수 분석 결과:")
    print("-" * 40)
    for i in range(n_modes):
        print(f"모드 {i+1}: {freq[i]:.3f} Hz ({omega[i]:.2f} rad/s)")

    return omega[:n_modes], eigenvectors[:, :n_modes]

# 3자유도 진동계
k1, k2, k3 = 1000, 1500, 800   # N/m
m1, m2, m3 = 1.0, 1.5, 0.8    # kg

K = np.array([
    [k1+k2,   -k2,    0],
    [-k2,   k2+k3,  -k3],
    [0,       -k3,   k3]
], dtype=float)

M = np.diag([m1, m2, m3])

omegas, modes = vibration_analysis(K, M)

# 모드 형상 시각화
positions = np.arange(1, 4)
fig, axes = plt.subplots(1, 3, figsize=(14, 5))

for i, ax in enumerate(axes):
    mode_shape = modes[:, i] / np.max(abs(modes[:, i]))
    ax.plot([0] + list(positions), [0] + list(mode_shape),
            'bo-', linewidth=2, markersize=10)
    ax.axhline(y=0, color='gray', linestyle='--', alpha=0.5)
    ax.set_title(f'모드 {i+1}\n{omegas[i]/(2*np.pi):.2f} Hz')
    ax.set_xlabel('자유도')
    ax.set_ylabel('정규화 모드 형상')
    ax.grid(True, alpha=0.3)
    ax.set_xticks([0, 1, 2, 3])
    ax.set_xticklabels(['고정단', '질량1', '질량2', '질량3'])

plt.suptitle('3자유도 진동계 모드 형상')
plt.tight_layout()
plt.savefig('vibration_modes.png', dpi=150)
plt.show()