import numpy as np
import matplotlib.pyplot as plt

# --- 段付き棒のパラメータ ---
# 3区間: [0, L1], [L1, L1+L2], [L1+L2, L1+L2+L3]
L = np.array([0.3, 0.4, 0.3])       # 各区間の長さ [m]
A = np.array([600e-6, 400e-6, 200e-6])  # 各区間の断面積 [m^2]
E = np.array([200e9, 200e9, 200e9])  # 各区間のヤング率 [Pa]

# 荷重: 区間1と2の境界に右向き80kN、区間2と3の境界に左向き30kN
P_ext = np.array([80e3, -30e3])  # 節点1, 節点2 に作用する外力 [N]

# --- 軸剛性の計算 ---
k = A * E / L  # 各区間の軸剛性 [N/m]

# --- 全体剛性行列の組み立て（変位の方法）---
# 内部節点: 節点1（区間1と2の境界）、節点2（区間2と3の境界）
# 全体剛性行列 K (2x2)
K = np.array([
    [k[0] + k[1], -k[1]],
    [-k[1], k[1] + k[2]]
])

# 連立方程式 K @ u = P_ext を解く
u = np.linalg.solve(K, P_ext)

print("=== 段付き棒の不静定解析 ===")
print(f"節点1の変位: {u[0]*1e6:.2f} μm")
print(f"節点2の変位: {u[1]*1e6:.2f} μm")

# --- 各区間の軸力と応力 ---
# 区間1: N1 = k[0] * u[0] （左壁から節点1への変位）
# 区間2: N2 = k[1] * (u[1] - u[0])
# 区間3: N3 = k[2] * (0 - u[1]) = -k[2] * u[1]
N = np.array([k[0]*u[0], k[1]*(u[1]-u[0]), -k[2]*u[1]])
sigma = N / A  # 各区間の応力 [Pa]

print(f"\n区間1: N = {N[0]/1e3:.2f} kN, σ = {sigma[0]/1e6:.2f} MPa")
print(f"区間2: N = {N[1]/1e3:.2f} kN, σ = {sigma[1]/1e6:.2f} MPa")
print(f"区間3: N = {N[2]/1e3:.2f} kN, σ = {sigma[2]/1e6:.2f} MPa")

# --- 反力 ---
R_A = -N[0]   # 左壁の反力（右向き正）
R_B = N[2]    # 右壁の反力（左向き正）
print(f"\n左壁の反力 R_A = {R_A/1e3:.2f} kN")
print(f"右壁の反力 R_B = {R_B/1e3:.2f} kN")
print(f"検算: R_A + 80 - 30 + R_B = {(R_A + 80e3 - 30e3 + R_B)/1e3:.4f} kN")

# --- 可視化 ---
fig, axes = plt.subplots(3, 1, figsize=(10, 10))
x_boundaries = np.cumsum(np.concatenate([[0], L]))

# (a) 応力分布
for i in range(3):
    axes[0].fill_between([x_boundaries[i], x_boundaries[i+1]],
                         sigma[i]/1e6, 0, alpha=0.3,
                         color='blue' if sigma[i] > 0 else 'red')
    axes[0].plot([x_boundaries[i], x_boundaries[i+1]],
                [sigma[i]/1e6, sigma[i]/1e6], 'k-', linewidth=2)
    axes[0].text((x_boundaries[i]+x_boundaries[i+1])/2, sigma[i]/1e6,
                f'{sigma[i]/1e6:.1f} MPa', ha='center', va='bottom' if sigma[i]>0 else 'top',
                fontsize=11, fontweight='bold')
axes[0].axhline(0, color='k', linewidth=0.5)
axes[0].set_ylabel('Stress [MPa]')
axes[0].set_title('(a) Stress Distribution')
axes[0].grid(True, alpha=0.3)

# (b) 軸力分布
for i in range(3):
    color = 'blue' if N[i] > 0 else 'red'
    axes[1].fill_between([x_boundaries[i], x_boundaries[i+1]],
                         N[i]/1e3, 0, alpha=0.3, color=color)
    axes[1].plot([x_boundaries[i], x_boundaries[i+1]],
                [N[i]/1e3, N[i]/1e3], 'k-', linewidth=2)
    axes[1].text((x_boundaries[i]+x_boundaries[i+1])/2, N[i]/1e3,
                f'{N[i]/1e3:.1f} kN', ha='center', va='bottom' if N[i]>0 else 'top',
                fontsize=11, fontweight='bold')
axes[1].axhline(0, color='k', linewidth=0.5)
axes[1].set_ylabel('Axial Force [kN]')
axes[1].set_title('(b) Axial Force Distribution')
axes[1].grid(True, alpha=0.3)

# (c) 変位分布
x_disp = [0]
u_disp = [0]
for i in range(3):
    x_end = x_boundaries[i+1]
    if i < 2:
        u_end = u[i]
    else:
        u_end = 0
    x_disp.append(x_end)
    u_disp.append(u_end)

axes[2].plot(x_disp, np.array(u_disp)*1e6, 'b-o', linewidth=2, markersize=8)
axes[2].axhline(0, color='k', linewidth=0.5)
axes[2].set_xlabel('Position x [m]')
axes[2].set_ylabel('Displacement [μm]')
axes[2].set_title('(c) Displacement Distribution')
axes[2].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('stepped_bar_indeterminate.png', dpi=150, bbox_inches='tight')
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# --- パラメータ ---
L1, L2 = 0.5, 0.5           # 区間の長さ [m]
A1, A2 = 500e-6, 300e-6     # 断面積 [m^2]
E1, E2 = 200e9, 70e9        # ヤング率 [Pa]（鋼 + アルミ）
alpha1, alpha2 = 12e-6, 23e-6  # 線膨張係数 [1/K]

# --- 温度変化の範囲 ---
dT = np.linspace(-100, 200, 300)

# --- 不静定応力の計算 ---
# N = -(alpha1*L1 + alpha2*L2)*dT / (L1/(A1*E1) + L2/(A2*E2))
numerator = -(alpha1 * L1 + alpha2 * L2) * dT
denominator = L1 / (A1 * E1) + L2 / (A2 * E2)
N = numerator / denominator

sigma1 = N / A1  # 区間1（鋼）の応力
sigma2 = N / A2  # 区間2（アルミ）の応力

# --- 変位（境界点の変位）---
u_boundary = alpha1 * dT * L1 + N * L1 / (A1 * E1)

# --- 可視化 ---
fig, axes = plt.subplots(1, 3, figsize=(16, 5))

# (1) 応力 vs 温度変化
axes[0].plot(dT, sigma1/1e6, 'b-', linewidth=2, label='Steel ($\\sigma_1$)')
axes[0].plot(dT, sigma2/1e6, 'r-', linewidth=2, label='Aluminum ($\\sigma_2$)')
axes[0].axhline(0, color='k', linewidth=0.5)
axes[0].axvline(0, color='k', linewidth=0.5)
axes[0].fill_between(dT, 0, sigma1/1e6, where=(sigma1<0), alpha=0.1, color='blue')
axes[0].fill_between(dT, 0, sigma1/1e6, where=(sigma1>0), alpha=0.1, color='blue')
axes[0].set_xlabel('Temperature Change $\\Delta T$ [K]')
axes[0].set_ylabel('Stress [MPa]')
axes[0].set_title('Thermal Stress in Stepped Bar')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# (2) 軸力 vs 温度変化
axes[1].plot(dT, N/1e3, 'k-', linewidth=2.5)
axes[1].axhline(0, color='k', linewidth=0.5)
axes[1].axvline(0, color='k', linewidth=0.5)
axes[1].fill_between(dT, 0, N/1e3, where=(N<0), alpha=0.2, color='red', label='Compression')
axes[1].fill_between(dT, 0, N/1e3, where=(N>0), alpha=0.2, color='blue', label='Tension')
axes[1].set_xlabel('Temperature Change $\\Delta T$ [K]')
axes[1].set_ylabel('Axial Force [kN]')
axes[1].set_title('Axial Force vs Temperature Change')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

# (3) 境界点変位 vs 温度変化
axes[2].plot(dT, u_boundary*1e6, 'g-', linewidth=2.5)
axes[2].axhline(0, color='k', linewidth=0.5)
axes[2].axvline(0, color='k', linewidth=0.5)
axes[2].set_xlabel('Temperature Change $\\Delta T$ [K]')
axes[2].set_ylabel('Boundary Displacement [μm]')
axes[2].set_title('Displacement at Section Boundary')
axes[2].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('thermal_indeterminate.png', dpi=150, bbox_inches='tight')
plt.show()

# --- 特定温度での数値確認 ---
dT_check = 100  # K
N_check = -(alpha1*L1 + alpha2*L2)*dT_check / denominator
print(f"ΔT = {dT_check} K のとき:")
print(f"  軸力 N = {N_check/1e3:.2f} kN")
print(f"  鋼の応力 σ1 = {N_check/A1/1e6:.2f} MPa")
print(f"  アルミの応力 σ2 = {N_check/A2/1e6:.2f} MPa")

import numpy as np
import matplotlib.pyplot as plt

# --- 両端固定梁 + 等分布荷重 ---
L = 4.0    # 梁の長さ [m]
w = 10e3   # 等分布荷重 [N/m]
E = 200e9  # ヤング率 [Pa]
I = 8.33e-6  # 断面二次モーメント [m^4]（H鋼200x100相当）
EI = E * I

# --- 解析解（両端固定梁の等分布荷重）---
# 反力: R_A = R_B = wL/2
# 固定端モーメント: M_A = M_B = -wL^2/12
# たわみ: v(x) = (w/(24EI)) * x^2 * (L-x)^2
R_A = w * L / 2
M_A = -w * L**2 / 12

x = np.linspace(0, L, 500)
v_analytical = (w / (24 * EI)) * x**2 * (L - x)**2

# --- 曲げモーメント分布 ---
M_analytical = M_A + R_A * x - w * x**2 / 2

# --- せん断力分布 ---
V_analytical = R_A - w * x

# --- 参考: 単純支持梁のたわみ ---
v_simple = (w / (24 * EI)) * x * (L**3 - 2*L*x**2 + x**3)
M_simple = w * L * x / 2 - w * x**2 / 2

# --- 可視化 ---
fig, axes = plt.subplots(3, 1, figsize=(10, 12))

# (a) たわみ曲線
axes[0].plot(x, v_analytical*1e3, 'b-', linewidth=2.5, label='Fixed-Fixed Beam')
axes[0].plot(x, v_simple*1e3, 'r--', linewidth=2, label='Simply Supported (reference)')
axes[0].set_ylabel('Deflection [mm]')
axes[0].set_title('(a) Deflection Curve')
axes[0].legend()
axes[0].grid(True, alpha=0.3)
axes[0].invert_yaxis()

# 最大たわみの注記
v_max_fixed = w * L**4 / (384 * EI)
v_max_simple = 5 * w * L**4 / (384 * EI)
axes[0].annotate(f'Max: {v_max_fixed*1e3:.3f} mm',
                xy=(L/2, v_max_fixed*1e3), xytext=(L/2+0.5, v_max_fixed*1e3+0.1),
                arrowprops=dict(arrowstyle='->', color='blue'), color='blue', fontsize=10)

# (b) 曲げモーメント図 (BMD)
axes[1].plot(x, M_analytical/1e3, 'b-', linewidth=2.5, label='Fixed-Fixed')
axes[1].plot(x, M_simple/1e3, 'r--', linewidth=2, label='Simply Supported')
axes[1].axhline(0, color='k', linewidth=0.5)
axes[1].fill_between(x, M_analytical/1e3, 0, alpha=0.15, color='blue')
axes[1].set_ylabel('Bending Moment [kN$\\cdot$m]')
axes[1].set_title('(b) Bending Moment Diagram (BMD)')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

# (c) せん断力図 (SFD)
axes[2].plot(x, V_analytical/1e3, 'b-', linewidth=2.5, label='Fixed-Fixed')
V_simple = w * L / 2 - w * x
axes[2].plot(x, V_simple/1e3, 'r--', linewidth=2, label='Simply Supported')
axes[2].axhline(0, color='k', linewidth=0.5)
axes[2].fill_between(x, V_analytical/1e3, 0, alpha=0.15, color='blue')
axes[2].set_xlabel('Position x [m]')
axes[2].set_ylabel('Shear Force [kN]')
axes[2].set_title('(c) Shear Force Diagram (SFD)')
axes[2].legend()
axes[2].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('fixed_beam_deflection.png', dpi=150, bbox_inches='tight')
plt.show()

print("=== 両端固定梁 vs 単純支持梁の比較 ===")
print(f"最大たわみ（両端固定）: {v_max_fixed*1e3:.4f} mm")
print(f"最大たわみ（単純支持）: {v_max_simple*1e3:.4f} mm")
print(f"たわみの比: {v_max_fixed/v_max_simple:.3f} （固定/単純 = 1/5）")
print(f"\n固定端モーメント: {M_A/1e3:.2f} kN·m")
print(f"中央モーメント（両端固定）: {(M_A + R_A*L/2 - w*(L/2)**2/2)/1e3:.2f} kN·m")
print(f"中央モーメント（単純支持）: {w*L**2/8/1e3:.2f} kN·m")

import numpy as np
import matplotlib.pyplot as plt

def solve_continuous_beam(spans, w_list, EI, n_elements_per_span=50):
    """
    連続梁の不静定解析（剛性法）

    Parameters
    ----------
    spans : list of float
        各スパンの長さ [m]
    w_list : list of float
        各スパンの等分布荷重 [N/m]
    EI : float
        曲げ剛性 [N·m^2]
    n_elements_per_span : int
        スパンあたりの要素数

    Returns
    -------
    x_global : array
        全体座標
    v_global : array
        たわみ分布
    M_global : array
        曲げモーメント分布
    reactions : array
        支持点反力
    """
    n_spans = len(spans)
    n_supports = n_spans + 1
    n_elem_total = n_elements_per_span * n_spans
    n_nodes = n_elem_total + 1

    # 全体長
    L_total = sum(spans)
    dx = L_total / n_elem_total

    # 支持点のノードインデックス
    support_nodes = [0]
    cumulative = 0
    for s in spans:
        cumulative += s
        node_idx = int(round(cumulative / dx))
        support_nodes.append(node_idx)

    # 全体剛性行列と荷重ベクトル（はり要素: Euler-Bernoulli）
    n_dof = 2 * n_nodes  # 各ノード: [v, θ]
    K_global = np.zeros((n_dof, n_dof))
    F_global = np.zeros(n_dof)

    for e in range(n_elem_total):
        Le = dx
        # 要素剛性行列
        ke = EI / Le**3 * np.array([
            [12,    6*Le,   -12,    6*Le],
            [6*Le,  4*Le**2, -6*Le, 2*Le**2],
            [-12,  -6*Le,    12,   -6*Le],
            [6*Le,  2*Le**2, -6*Le, 4*Le**2]
        ])

        # 等分布荷重の等価節点力
        x_mid = (e + 0.5) * dx
        # どのスパンに属するか
        cum = 0
        w_e = 0
        for s_idx, s_len in enumerate(spans):
            if x_mid < cum + s_len + 1e-10:
                w_e = w_list[s_idx]
                break
            cum += s_len

        fe = w_e * Le / 2 * np.array([1, Le/6, 1, -Le/6])

        # 自由度マッピング
        dofs = [2*e, 2*e+1, 2*(e+1), 2*(e+1)+1]
        for i in range(4):
            F_global[dofs[i]] += fe[i]
            for j in range(4):
                K_global[dofs[i], dofs[j]] += ke[i, j]

    # 境界条件: 支持点でv=0（ピン支持 / ローラー支持）
    # 拘束自由度: 各支持点のv（鉛直変位）
    constrained_dofs = [2 * node for node in support_nodes]

    free_dofs = [i for i in range(n_dof) if i not in constrained_dofs]

    K_ff = K_global[np.ix_(free_dofs, free_dofs)]
    F_f = F_global[free_dofs]

    u_f = np.linalg.solve(K_ff, F_f)

    # 全変位ベクトルの組み立て
    u_global = np.zeros(n_dof)
    for i, dof in enumerate(free_dofs):
        u_global[dof] = u_f[i]

    # たわみの抽出
    x_global = np.linspace(0, L_total, n_nodes)
    v_global = u_global[0::2]  # 偶数インデックス = たわみ

    # 反力の計算
    F_reaction = K_global @ u_global - F_global
    reactions = np.array([F_reaction[2*node] for node in support_nodes])

    # 曲げモーメントの計算（差分近似）
    theta = u_global[1::2]  # 奇数インデックス = 回転角
    M_global = np.zeros(n_nodes)
    for e in range(n_elem_total):
        Le = dx
        # 要素端モーメント
        u_e = np.array([v_global[e], theta[e], v_global[e+1], theta[e+1]])
        ke = EI / Le**3 * np.array([
            [12,    6*Le,   -12,    6*Le],
            [6*Le,  4*Le**2, -6*Le, 2*Le**2],
            [-12,  -6*Le,    12,   -6*Le],
            [6*Le,  2*Le**2, -6*Le, 4*Le**2]
        ])
        f_e = ke @ u_e
        # 要素左端と右端のモーメント
        M_global[e] += f_e[1] / 2 if e > 0 else f_e[1]
        M_global[e+1] += f_e[3] / 2 if e < n_elem_total - 1 else f_e[3]

    return x_global, v_global, M_global, reactions, support_nodes

# --- 2スパン連続梁 ---
EI = 200e9 * 8.33e-6  # [N·m^2]
spans_2 = [5.0, 5.0]  # 等スパン
w_2 = [10e3, 10e3]     # 等分布荷重 [N/m]

x2, v2, M2, R2, sn2 = solve_continuous_beam(spans_2, w_2, EI)

# --- 3スパン連続梁 ---
spans_3 = [4.0, 5.0, 4.0]
w_3 = [10e3, 10e3, 10e3]

x3, v3, M3, R3, sn3 = solve_continuous_beam(spans_3, w_3, EI)

# --- 可視化 ---
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 2スパン: たわみ
axes[0,0].plot(x2, v2*1e3, 'b-', linewidth=2)
for node in sn2:
    axes[0,0].plot(x2[node], 0, 'k^', markersize=12)
axes[0,0].axhline(0, color='k', linewidth=0.5)
axes[0,0].set_title('2-Span Continuous Beam: Deflection')
axes[0,0].set_ylabel('Deflection [mm]')
axes[0,0].invert_yaxis()
axes[0,0].grid(True, alpha=0.3)

# 2スパン: BMD
axes[0,1].plot(x2, M2/1e3, 'b-', linewidth=2)
axes[0,1].fill_between(x2, M2/1e3, 0, alpha=0.15, color='blue')
axes[0,1].axhline(0, color='k', linewidth=0.5)
axes[0,1].set_title('2-Span Continuous Beam: BMD')
axes[0,1].set_ylabel('Bending Moment [kN$\\cdot$m]')
axes[0,1].grid(True, alpha=0.3)

# 3スパン: たわみ
axes[1,0].plot(x3, v3*1e3, 'r-', linewidth=2)
for node in sn3:
    axes[1,0].plot(x3[node], 0, 'k^', markersize=12)
axes[1,0].axhline(0, color='k', linewidth=0.5)
axes[1,0].set_title('3-Span Continuous Beam: Deflection')
axes[1,0].set_xlabel('Position x [m]')
axes[1,0].set_ylabel('Deflection [mm]')
axes[1,0].invert_yaxis()
axes[1,0].grid(True, alpha=0.3)

# 3スパン: BMD
axes[1,1].plot(x3, M3/1e3, 'r-', linewidth=2)
axes[1,1].fill_between(x3, M3/1e3, 0, alpha=0.15, color='red')
axes[1,1].axhline(0, color='k', linewidth=0.5)
axes[1,1].set_title('3-Span Continuous Beam: BMD')
axes[1,1].set_xlabel('Position x [m]')
axes[1,1].set_ylabel('Bending Moment [kN$\\cdot$m]')
axes[1,1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('continuous_beam.png', dpi=150, bbox_inches='tight')
plt.show()

# --- 反力の表示 ---
print("=== 2スパン連続梁の反力 ===")
support_labels = ['A', 'B', 'C']
for i, r in enumerate(R2):
    print(f"  R_{support_labels[i]} = {r/1e3:.2f} kN")
print(f"  合計: {sum(R2)/1e3:.2f} kN (外力合計: {sum(w_2[j]*spans_2[j] for j in range(len(spans_2)))/1e3:.2f} kN)")

print("\n=== 3スパン連続梁の反力 ===")
support_labels_3 = ['A', 'B', 'C', 'D']
for i, r in enumerate(R3):
    print(f"  R_{support_labels_3[i]} = {r/1e3:.2f} kN")
print(f"  合計: {sum(R3)/1e3:.2f} kN (外力合計: {sum(w_3[j]*spans_3[j] for j in range(len(spans_3)))/1e3:.2f} kN)")

import numpy as np
import matplotlib.pyplot as plt

class BeamSolver:
    """Euler-Bernoulliはりの不静定解析ソルバー（剛性法）"""

    def __init__(self, L, EI, n_elem=100):
        """
        Parameters
        ----------
        L : float  梁の全長 [m]
        EI : float 曲げ剛性 [N·m^2]
        n_elem : int 要素数
        """
        self.L = L
        self.EI = EI
        self.n_elem = n_elem
        self.n_nodes = n_elem + 1
        self.dx = L / n_elem
        self.n_dof = 2 * self.n_nodes
        self.supports = {}      # {node: 'pin', 'roller', 'fixed'}
        self.point_loads = []   # [(node, force)]
        self.dist_loads = []    # [(x_start, x_end, w)]

    def add_support(self, x, stype='pin'):
        """支持条件の追加"""
        node = int(round(x / self.dx))
        self.supports[node] = stype

    def add_point_load(self, x, P):
        """集中荷重の追加"""
        node = int(round(x / self.dx))
        self.point_loads.append((node, P))

    def add_distributed_load(self, x_start, x_end, w):
        """等分布荷重の追加"""
        self.dist_loads.append((x_start, x_end, w))

    def solve(self):
        """剛性法による求解"""
        dx = self.dx
        K = np.zeros((self.n_dof, self.n_dof))
        F = np.zeros(self.n_dof)

        # 全体剛性行列の組み立て
        for e in range(self.n_elem):
            ke = self.EI / dx**3 * np.array([
                [12,    6*dx,   -12,    6*dx],
                [6*dx,  4*dx**2, -6*dx, 2*dx**2],
                [-12,  -6*dx,    12,   -6*dx],
                [6*dx,  2*dx**2, -6*dx, 4*dx**2]
            ])
            dofs = [2*e, 2*e+1, 2*(e+1), 2*(e+1)+1]
            for i in range(4):
                for j in range(4):
                    K[dofs[i], dofs[j]] += ke[i, j]

        # 集中荷重
        for node, P in self.point_loads:
            F[2*node] += P

        # 等分布荷重
        for x_s, x_e, w in self.dist_loads:
            for e in range(self.n_elem):
                x_mid = (e + 0.5) * dx
                if x_s <= x_mid <= x_e:
                    fe = w * dx / 2 * np.array([1, dx/6, 1, -dx/6])
                    dofs = [2*e, 2*e+1, 2*(e+1), 2*(e+1)+1]
                    for i in range(4):
                        F[dofs[i]] += fe[i]

        # 境界条件の適用
        constrained = []
        for node, stype in self.supports.items():
            constrained.append(2*node)  # たわみ = 0
            if stype == 'fixed':
                constrained.append(2*node + 1)  # 回転角 = 0

        free = [i for i in range(self.n_dof) if i not in constrained]
        u = np.zeros(self.n_dof)
        u[free] = np.linalg.solve(K[np.ix_(free, free)], F[free])

        # 反力
        reactions = K @ u - F

        # 結果の格納
        self.x = np.linspace(0, self.L, self.n_nodes)
        self.v = u[0::2]
        self.theta = u[1::2]
        self.reactions = {node: reactions[2*node] for node in self.supports}
        self.moment_reactions = {node: reactions[2*node+1]
                                 for node, s in self.supports.items() if s == 'fixed'}

        # 曲げモーメントの計算
        self.M = np.zeros(self.n_nodes)
        for e in range(self.n_elem):
            ue = np.array([self.v[e], self.theta[e], self.v[e+1], self.theta[e+1]])
            ke = self.EI / dx**3 * np.array([
                [12,    6*dx,   -12,    6*dx],
                [6*dx,  4*dx**2, -6*dx, 2*dx**2],
                [-12,  -6*dx,    12,   -6*dx],
                [6*dx,  2*dx**2, -6*dx, 4*dx**2]
            ])
            fe = ke @ ue
            if e == 0:
                self.M[e] = fe[1]
            else:
                self.M[e] = (self.M[e] + fe[1]) / 2 if self.M[e] != 0 else fe[1]
            self.M[e+1] = fe[3]

        return self

# --- 応用例: プロッツ梁（一端固定・他端ピン + 集中荷重）---
L = 6.0
EI = 200e9 * 8.33e-6

beam = BeamSolver(L, EI, n_elem=200)
beam.add_support(0, 'fixed')
beam.add_support(L, 'pin')
beam.add_point_load(L/3, -50e3)   # 下向き50kN

result = beam.solve()

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

axes[0].plot(result.x, result.v*1e3, 'b-', linewidth=2.5)
axes[0].axhline(0, color='k', linewidth=0.5)
axes[0].set_ylabel('Deflection [mm]')
axes[0].set_title('Propped Cantilever: Deflection')
axes[0].invert_yaxis()
axes[0].grid(True, alpha=0.3)

axes[1].plot(result.x, result.M/1e3, 'b-', linewidth=2.5)
axes[1].fill_between(result.x, result.M/1e3, 0, alpha=0.15, color='blue')
axes[1].axhline(0, color='k', linewidth=0.5)
axes[1].set_xlabel('Position x [m]')
axes[1].set_ylabel('Bending Moment [kN$\\cdot$m]')
axes[1].set_title('Propped Cantilever: BMD')
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('propped_cantilever.png', dpi=150, bbox_inches='tight')
plt.show()

print("=== プロッツ梁の反力 ===")
for node, r in result.reactions.items():
    x_pos = node * result.dx
    print(f"  x = {x_pos:.1f} m: R = {r/1e3:.2f} kN")
for node, m in result.moment_reactions.items():
    x_pos = node * result.dx
    print(f"  x = {x_pos:.1f} m: M = {m/1e3:.2f} kN·m")