d0 = [0,0];
v0 = [2,3];
a0 = [0,0];
dt = 1e-3;
n = 9001;
Mass = @(d,v,a) mass(d,v,a);
Damping = @(d,v,a) damping(d,v,a);
Stiffness = @(d,v,a) stiffness(d,v,a);
Force = @(t) force(t);
[dM,vM,aM] = Galpha_integration(dt,n,d0,v0,a0,Mass,Damping,Stiffness,Force);
tspan = 0:dt:(n-1)*dt;
plot(tspan,dM(:,1),'b','linewidth',2.5);
hold on
plot(tspan,sin(2*tspan),'r','linewidth',1.5) ;
xlabel('t');
ylabel('u_1')
legend('Generalized-alpha','Exact')
figure(2)
plot(tspan,dM(:,2),'b','linewidth',2.5);
hold on
plot(tspan,sin(3*tspan),'r','linewidth',1.5)
legend('Generalized-alpha','Exact')
xlabel('t');
ylabel('u_2')
function [dM,vM,aM] = Galpha_integration(dt,n,d0,v0,a0,Mass,Damp,Stiffness,Force)
alpha_m = 1/2; alpha_f = 1/2; gamma = 1/2-alpha_m+alpha_f; beta = 1/4*(1-alpha_m+alpha_f)^2;
d = length(d0); dM = zeros(n,d); vM = zeros(n,d); aM = zeros(n,d); dM(1,:) = d0; vM(1,:) = v0; aM(1,:) = a0;
options = optimoptions('fsolve','Algorithm','levenberg-marquardt','Display','None');
for i = 1:n-1
dc = dM(i,:); vc = vM(i,:); ac = aM(i,:); tc = (i-1)*dt;
fun = @(x) Generalized_alpha(x,dc,vc,ac,dt,tc,Mass,Damp,Stiffness,Force,alpha_m,alpha_f);
y = fsolve(fun,ac,options); y = (y(:))'; dn = dc+dt*vc+dt^2/2*((1-2*beta)*ac+2*beta*y);
vn = vc+dt*((1-gamma)*ac+gamma*y); dM(i+1,:) = dn; vM(i+1,:) = vn; aM(i+1,:) = y;
end
end
function equ = Generalized_alpha(x,dc,vc,ac,dt,tc,Mass,Damp,Stiffness,Force,alpha_m,alpha_f)
gamma = 1/2-alpha_m+alpha_f; beta = 1/4*(1-alpha_m+alpha_f)^2;
dn = dc+dt*vc+dt^2/2*((1-2*beta)*ac+2*beta*x); vn = vc+dt*((1-gamma)*ac+gamma*x);
df = (1-alpha_f)*dn+alpha_f*dc; vf = (1-alpha_f)*vn+alpha_f*vc; am = (1-alpha_m)*x+alpha_m*ac;
tf = (1-alpha_f)*(tc+dt)+alpha_f*tc; M = Mass(df,vf,am); C = Damp(df,vf,am); K = Stiffness(df,vf,am);
Ff = Force(tf); equ = M*am'+C*vf'+K*df'-(Ff(:))';
end
function M = mass(d,v,a)
M = [1,0;0,1];
end
function C = damping(d,v,a)
C = zeros(2,2);
end
function K = stiffness(d,v,a)
K = [4,0;0,9];
end
function F = force(t)
F = zeros(2,1);
end
