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);

gamma = 0.5; beta = 0.25;   % ortalama ivme (kosulsuz kararli Newmark semasi)

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);   % serbest titresim: dis kuvvet yok

    % Newmark-beta integrasyonu - dogrudan satir ici (fonksiyon yok, cunku
    % betik-ici yerel fonksiyonlar Octave'da .m dosyasi olarak DOGRUDAN
    % calistirildiginda calismiyor - bu platform tarafindan Octave ile
    % test edilip dogrulandi, 05.08.2026)
    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

    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