Ana Sayfa / MATLAB / Python / Problem 5

Newmark-β Zaman Tanım Alanı İntegrasyonu

Sismik analizin sayısal çekirdeği. Burada serbest titreşimle sınıyoruz — çünkü serbest titreşimin tam analitik çözümü var, yani integratörün hatasını mutlak olarak ölçebiliyoruz.

Neden serbest titreşimle test ediyoruz? Gerçek bir deprem kaydıyla test etseydik, "doğru cevabı" bilemezdik — sadece iki sayısal sonucu karşılaştırabilirdik. Serbest titreşimde ise analitik çözüm x(t)=x₀e−ζωnt(cos ωdt + (ζωn/ωd)sin ωdt) tam olarak bilinir. Bir integratörü doğrulamanın tek dürüst yolu budur.

PYTHON newmark_beta.py

import numpy as np

def newmark_sdof(m, c, k, p, dt, u0=0, v0=0, gamma=0.5, beta=0.25):
    """Newmark-beta (varsayilan: ortalama ivme, kosulsuz kararli)"""
    n = len(p)
    u = np.zeros(n); v = np.zeros(n); a = np.zeros(n)
    u[0], v[0] = u0, v0
    a[0] = (p[0] - c*v0 - k*u0) / m

    a1 = m/(beta*dt**2) + gamma*c/(beta*dt)
    a2 = m/(beta*dt) + (gamma/beta - 1)*c
    a3 = (1/(2*beta) - 1)*m + dt*(gamma/(2*beta) - 1)*c
    k_hat = k + a1

    for i in range(n-1):
        p_hat = p[i+1] + a1*u[i] + a2*v[i] + a3*a[i]
        u[i+1] = p_hat / k_hat
        v[i+1] = (gamma/(beta*dt))*(u[i+1]-u[i]) + (1-gamma/beta)*v[i] \
                 + dt*(1-gamma/(2*beta))*a[i]
        a[i+1] = (u[i+1]-u[i])/(beta*dt**2) - v[i]/(beta*dt) \
                 - (1/(2*beta)-1)*a[i]
    return u, v, a

# Bilinen sistem: analitik cozumu TAM olarak bilinen serbest titresim
m, wn, zeta = 1000.0, 10.0, 0.05
k = m*wn**2
c = 2*zeta*wn*m
u0, v0 = 0.01, 0.0
wd = wn*np.sqrt(1 - zeta**2)

print("dt(s)    T/dt    maks hata(%)   hata orani")
prev_err = None
for dt in [0.05, 0.02, 0.01, 0.005, 0.002]:
    T_end = 3.0
    n = int(T_end/dt) + 1
    t = np.linspace(0, T_end, n)
    p = np.zeros(n)                       # serbest titresim: dis kuvvet yok

    u, _, _ = newmark_sdof(m, c, k, p, dt, u0, v0)

    # Analitik cozum
    ua = u0*np.exp(-zeta*wn*t)*(np.cos(wd*t) + (zeta*wn/wd)*np.sin(wd*t))
    err = np.max(np.abs(u - ua))/u0*100

    ratio = f"{prev_err/err:.2f}x" if prev_err else "-"
    print(f"{dt:.3f}   {2*np.pi/wn/dt:6.1f}   {err:9.4f}     {ratio}")
    prev_err = err
dt(s) T/dt maks hata(%) hata orani 0.050 12.6 15.1198 - 0.020 31.4 2.4519 6.17x 0.010 62.8 0.6130 4.00x 0.005 125.7 0.1534 3.99x 0.002 314.2 0.0246 6.25x
✓ Bu platform tarafından gerçekten çalıştırıldı — analitik çözümle karşılaştırıldı
Sayısal Bulgu: İkinci Mertebeden Yakınsama Zaman adımı yarıya indiğinde hata dörde bölünüyor (0.020→0.010 geçişinde tam 4.00×, 0.010→0.005'te 3.99×). Bu, Newmark-β ortalama ivme yönteminin O(Δt²) — yani ikinci mertebeden doğru — olduğunun sayısal kanıtıdır. Bu deseni kendi kodunuzda göremiyorsanız, ya integratörünüzde bir hata vardır ya da yanlışlıkla birinci mertebeden bir şema kullanıyorsunuzdur.

Pratik kural: T/Δt ≥ 20 (yani periyot başına en az 20 adım) genellikle yeterlidir — tabloda T/Δt=31.4'te hata zaten %2.5, T/Δt=62.8'te %0.6.

MATLAB newmark_beta.m

m = 1000.0; wn = 10.0; zeta = 0.05;
k = m*wn^2;
c = 2*zeta*wn*m;
u0 = 0.01; v0 = 0.0;
wd = wn*sqrt(1 - zeta^2);

fprintf('dt(s)    T/dt    maks hata(%%)\n');
for dt = [0.05, 0.02, 0.01, 0.005, 0.002]
    T_end = 3.0;
    n = round(T_end/dt) + 1;
    t = linspace(0, T_end, n);
    p = zeros(1, n);

    u = newmark_sdof(m, c, k, p, dt, u0, v0);

    ua = u0*exp(-zeta*wn*t).*(cos(wd*t) + (zeta*wn/wd)*sin(wd*t));
    err = max(abs(u - ua))/u0*100;

    fprintf('%.3f   %6.1f   %9.4f\n', dt, 2*pi/wn/dt, err);
end

% --- Yerel fonksiyon (MATLAB kurali: dosyanin EN SONUNDA) ---
function u = newmark_sdof(m, c, k, p, dt, u0, v0)
    gamma = 0.5; beta = 0.25;   % ortalama ivme
    n = length(p);
    u = zeros(1,n); v = zeros(1,n); a = zeros(1,n);
    u(1) = u0; v(1) = v0;
    a(1) = (p(1) - c*v0 - k*u0)/m;

    a1 = m/(beta*dt^2) + gamma*c/(beta*dt);
    a2 = m/(beta*dt) + (gamma/beta - 1)*c;
    a3 = (1/(2*beta) - 1)*m + dt*(gamma/(2*beta) - 1)*c;
    k_hat = k + a1;

    for i = 1:n-1
        p_hat = p(i+1) + a1*u(i) + a2*v(i) + a3*a(i);
        u(i+1) = p_hat/k_hat;
        v(i+1) = (gamma/(beta*dt))*(u(i+1)-u(i)) + (1-gamma/beta)*v(i) ...
                 + dt*(1-gamma/(2*beta))*a(i);
        a(i+1) = (u(i+1)-u(i))/(beta*dt^2) - v(i)/(beta*dt) ...
                 - (1/(2*beta)-1)*a(i);
    end
end
dt(s) T/dt maks hata(%) 0.050 12.6 15.1198 0.020 31.4 2.4519 0.010 62.8 0.6130 0.005 125.7 0.1534 0.002 314.2 0.0246

⚠ Bu kod gerçek bir MATLAB'da çalıştırılmadı. Python'daki aynı algoritma, MATLAB'ın 1-indeksli dizileri ve eleman-bazlı çarpım operatörü (.*) kuralına dikkat edilerek yazıldı. Yerel fonksiyon, MATLAB kuralı gereği dosyanın en sonuna konuldu. Kendi MATLAB'ınızda çalıştırıp doğrulayın.

← Problem 4 SHM: Aynı Sinyalden Frekans Çıkarımı →