From f6704e3029fd22254e6d8fa4ab33ba51f3985b67 Mon Sep 17 00:00:00 2001 From: jain-shobhit <22640963+jain-shobhit@users.noreply.github.com> Date: Tue, 8 Sep 2020 20:39:09 +0200 Subject: [PATCH] Modeselection (#2) * Modeselection * minor changes #1 * ensuring sparse/full eigenvalue computation in the reduced/full case. * Commentary Co-authored-by: gergelybuza <57031910+gergelybuza@users.noreply.github.com> --- @SSR/DF_P.m | 3 +- @SSR/FRC.m | 12 ++ @SSR/SSM2.m | 53 +++++++ @SSR/SSMcomps.m | 15 ++ @SSR/SSR.m | 95 ++++++++++-- examples/ModifiedShawPierrePeriodic.m | 7 + examples/ModifiedShawPierreQuasiPeriodic.m | 7 + examples/motiv_example.m | 101 +++++++++++++ .../ndof_oscillator_chain/OscillatorChain.m | 6 + examples/nonsmooth/NonSmoothOscillator.m | 7 + examples/vonKarman/example_curvedbeam.m | 99 +++++++++++++ examples/vonKarman/example_straightbeam.m | 136 ++++++++++++++++++ .../vonKarman/model/AssembleLinearStiffness.m | 33 +++++ examples/vonKarman/model/AssembleMass.m | 30 ++++ .../model/AssembleTangentStiffness.m | 35 +++++ examples/vonKarman/model/Beam_Model.m | 89 ++++++++++++ examples/vonKarman/model/BoundaryCoundition.m | 11 ++ examples/vonKarman/model/ElementMass.m | 22 +++ examples/vonKarman/model/ElementStiffness.m | 34 +++++ examples/vonKarman/model/Geometry.m | 9 ++ .../vonKarman/model/NonlinearityJacobianVK.m | 16 +++ examples/vonKarman/model/NonlinearityVK.m | 14 ++ examples/vonKarman/model/PlotDeformation.m | 25 ++++ examples/vonKarman/model/S11.m | 80 +++++++++++ .../vonKarman/model/compute_rotation_matrix.m | 13 ++ examples/vonKarman/model/sortVMsVK.m | 11 ++ modeselection/S11reduct.m | 16 +++ modeselection/copyfromfull.m | 15 ++ modeselection/modeselect.m | 101 +++++++++++++ modeselection/scalg.m | 24 ++++ 30 files changed, 1104 insertions(+), 15 deletions(-) create mode 100644 @SSR/FRC.m create mode 100644 @SSR/SSM2.m create mode 100644 @SSR/SSMcomps.m create mode 100644 examples/motiv_example.m create mode 100644 examples/vonKarman/example_curvedbeam.m create mode 100644 examples/vonKarman/example_straightbeam.m create mode 100644 examples/vonKarman/model/AssembleLinearStiffness.m create mode 100644 examples/vonKarman/model/AssembleMass.m create mode 100644 examples/vonKarman/model/AssembleTangentStiffness.m create mode 100644 examples/vonKarman/model/Beam_Model.m create mode 100644 examples/vonKarman/model/BoundaryCoundition.m create mode 100644 examples/vonKarman/model/ElementMass.m create mode 100644 examples/vonKarman/model/ElementStiffness.m create mode 100644 examples/vonKarman/model/Geometry.m create mode 100644 examples/vonKarman/model/NonlinearityJacobianVK.m create mode 100644 examples/vonKarman/model/NonlinearityVK.m create mode 100644 examples/vonKarman/model/PlotDeformation.m create mode 100644 examples/vonKarman/model/S11.m create mode 100644 examples/vonKarman/model/compute_rotation_matrix.m create mode 100644 examples/vonKarman/model/sortVMsVK.m create mode 100644 modeselection/S11reduct.m create mode 100644 modeselection/copyfromfull.m create mode 100644 modeselection/modeselect.m create mode 100644 modeselection/scalg.m diff --git a/@SSR/DF_P.m b/@SSR/DF_P.m index abca5bd..e844e13 100644 --- a/@SSR/DF_P.m +++ b/@SSR/DF_P.m @@ -14,5 +14,6 @@ D = blkdiag(DSC{:}); BU = blkdiag(BU{:}); end -DF_P = speye( O.nt * O.n ) + BU * get_ConvMtx(O) * D; +DD = BU * get_ConvMtx(O) * D; +DF_P = speye( size(DD) ) + DD; end \ No newline at end of file diff --git a/@SSR/FRC.m b/@SSR/FRC.m new file mode 100644 index 0000000..0434af1 --- /dev/null +++ b/@SSR/FRC.m @@ -0,0 +1,12 @@ +function [om, N] = FRC(O,T_range,ncontsteps) +% function that computes the FRC of S for a given range + +[om,cont_sol,~] = O.sequential_continuation(T_range,'ncontsteps',ncontsteps); + +N = nan(length(om),1); +for j = 1:length(om) + if ~isempty(cont_sol{j}) + N(j) = O.dt*norm(cont_sol{j},'fro'); + end +end +end \ No newline at end of file diff --git a/@SSR/SSM2.m b/@SSR/SSM2.m new file mode 100644 index 0000000..be074bd --- /dev/null +++ b/@SSR/SSM2.m @@ -0,0 +1,53 @@ +function W = SSM2(O,S) +% A faster computation of the leading order term of the SSM + +n = length(O.M); +m = O.n; +W = cell(1,n-m); + +dim = (2*m)^2; + +A = [zeros(m),eye(m);-O.K1,-O.C1]; % linear coefficient of master system +B0 = zeros(dim); + +delta = eye(2*m); % kronecker delta + +% assembly of the 2D equivalent to tensor B as a part of vectorizing the +% problem (first half that is independent of enslaved mode) +for i = 1:dim + for j = 1:dim + q = floor((j-1)/2/m)+1; + t = floor((i-1)/2/m)+1; + r = j - (q-1) * 2 * m; + s = i - (t-1) * 2 * m; + B0(i,j) = 2 * A(r,s) * A(q,t) +... + A(q,:) * A(:,t) * delta(r,s) +... + A(r,:) * A(:,s) * delta(q,t); % eq. (2.14) in thesis + end +end + + +% Computing SSM coefficient Wk for each enslaved mode +for k = 1:n-m + B1 = zeros(dim); + Sk = zeros(2*m); + Sk(1:m,1:m) = S{k}; % this is R_k, nonlinear quadratic coeffs in original eq. + + % Computing second half of B + for i = 1:dim + for j = 1:dim + q = floor((j-1)/2/m)+1; + t = floor((i-1)/2/m)+1; + r = j - (q-1) * 2 * m; + s = i - (t-1) * 2 * m; + B1(i,j) = 2 * (O.zeta2(k)) * (O.omega2(k)) *... + (A(q,t) * delta(r,s) + A(r,s) * delta(q,t)); + end + end + + B = B0 + B1 + (O.omega2(k))^2 * eye(dim); + Wk = -B\Sk(:); % quadratic SSM coeffs in vector form + W{k} = reshape(Wk,size(Sk)); % converting SSM coeffs back to matrix form +end + +end \ No newline at end of file diff --git a/@SSR/SSMcomps.m b/@SSR/SSMcomps.m new file mode 100644 index 0000000..ea4b408 --- /dev/null +++ b/@SSR/SSMcomps.m @@ -0,0 +1,15 @@ +function SSMcomps(O) +% computes additional material required for SSM computation +[V1, dd] = eig(full(O.K),full(O.M),'vector'); +[~, ind] = sort(dd); +V1 = V1(:,ind); +V2 = V1(:,setdiff(1:length(O.M),O.mode_choice)); +mu2 = diag(V2.' * O.M * V2); +U2 = V2 * diag( 1./ sqrt(mu2) ); +omega2 = sqrt(diag((U2.' * O.K * U2))); +zeta2 = diag((U2.' * O.C * U2))./ (2*omega2); +O.U2store = U2; +O.omega2store = omega2; +O.zeta2store = zeta2; +O.SSMdone = 1; +end \ No newline at end of file diff --git a/@SSR/SSR.m b/@SSR/SSR.m index 362b5e4..ac9c84f 100644 --- a/@SSR/SSR.m +++ b/@SSR/SSR.m @@ -1,6 +1,15 @@ classdef SSR < handle - %SSR Summary of this class goes here - % Detailed explanation goes here + %SSR (SteadyStateResponse) class implements the integral equations + %approach to the computation of the steady-state response in + %(quasi)periodically-forced dynamical systems. + + % The theory of computation (without model reduction) is decribed in + % the following article: + % + + % The theory for model reduction is described in the following article: + % G. Buza, S. Jain, G. Haller, Using Spectral Submanifolds for Optimal + % Mode Selection in Model Reduction, (2020) Preprint available on arXiv.org properties M % Mass matrix @@ -35,6 +44,8 @@ DLvec % Derivative dLdT Jvec % Green's function for velocities (calculated over -T:dt:T) ConvMtx % convolution matrix for displacements + A % reformulated convolution matrix + Lint % precomputed integrals of Green's functions isupdated % check if different quantitities are updated according to T weights % Weights vector for higher order integration nt % number of nodes in the time mesh @@ -51,30 +62,80 @@ Evinv % The inverse Fourier transform exponential matrix Ekron % The Fourier transform exponential matrix kronecker product with Identity Einvkron % The inverseFourier transform exponential matrix kronecker product with Identity + mode_choice % Selection of modes specified by the user + + % required only for the SSM computation in the mode selection + % procedure + zeta % The damping coefficients of the master modes + omega2 % The eigenfrequencies of the enslaved modes + zeta2 % The damping coefficients of the enslaved modes + K1 % Stiffness matrix of the enslaved modes (in modal coordinates) + C1 % Damping matrix of the enslaved modes (in modal coordinates) + U2 % Enslaved eigenmodes + % storage for some of these + U2store + zeta2store + omega2store + SSMdone % Boolean that ensures that the eigenvalues are only computed once end methods - function O = SSR(M,C,K) + function O = SSR(M,C,K,varargin) % Basic system properties O.M = M; O.C = C; - O.K = K; - O.n = length(M); + O.K = K; + if nargin == 3 + n = length(M); + O.mode_choice = 1:n; + [VV, dd] = eig(full(K),full(M)); + else + O.mode_choice = varargin{1}; + [VV, dd] = eigs(K,M,max(O.mode_choice),'SM'); + end + O.n = length(O.mode_choice); - [V, ~] = eig(full(K),full(M)); + dd = diag(dd); + [~, ind] = sort(dd); + VV = VV(:,ind); + V = VV(:,O.mode_choice); mu = diag(V.' * M * V); O.U = V * diag( 1./ sqrt(mu) ); % mass normalized VMs O.omega0 = sqrt(diag((O.U.' * K * O.U))); - zeta = diag((O.U.' * C * O.U))./ (2*O.omega0); - O.o = find(zeta>1); - O.c = find(zeta==1); - O.u = find(zeta<1); - lambda1 = (-zeta + sqrt(zeta.^2 - 1)).* O.omega0; - lambda2 = (-zeta - sqrt(zeta.^2 - 1)).* O.omega0; + O.K1 = O.U.' * K * O.U; + O.zeta = diag((O.U.' * C * O.U))./ (2*O.omega0); + O.C1 = O.U.' * C * O.U; + O.o = find(O.zeta>1); + O.c = find(O.zeta==1); + O.u = find(O.zeta<1); + lambda1 = (-O.zeta + sqrt(O.zeta.^2 - 1)).* O.omega0; + lambda2 = (-O.zeta - sqrt(O.zeta.^2 - 1)).* O.omega0; O.alpha = real(lambda2); O.omega = abs(imag(lambda2)); O.beta = lambda1; O.gamma = lambda2; + O.SSMdone = 0; + end + + function U2 = get.U2(O) + if ~O.SSMdone + SSMcomps(O); + end + U2 = O.U2store; + end + + function omega2 = get.omega2(O) + if ~O.SSMdone + SSMcomps(O); + end + omega2 = O.omega2store; + end + + function zeta2 = get.zeta2(O) + if ~O.SSMdone + SSMcomps(O); + end + zeta2 = O.zeta2store; end function set.T(O,T) @@ -188,7 +249,7 @@ Jvec = O.Jvec; end end - + update_Jvec(O) J = J(O,t,T) @@ -224,7 +285,7 @@ DF_Q = DF_Q(O,x) - [x,xd] = LinearResponse(O) + [x, xd] = LinearResponse(O) [x, xd] = Picard(O,varargin) @@ -233,6 +294,12 @@ [x0,tol,maxiter] = parse_iteration_inputs(O,inputs) [OMEGA, SOL, PICARD] = sequential_continuation(O,range,varargin) + + W = SSM2(O,S) + + [U2,omega2,zeta2] = SSMcomps(O) + + [om, N] = FRC(O,T_range,ncontsteps) end methods(Static) diff --git a/examples/ModifiedShawPierrePeriodic.m b/examples/ModifiedShawPierrePeriodic.m index d3c38c1..b4fa3f7 100644 --- a/examples/ModifiedShawPierrePeriodic.m +++ b/examples/ModifiedShawPierrePeriodic.m @@ -2,6 +2,13 @@ % m*x1dd + k(x1-x2) + c(x1d-x2d) + g*x1^3 = f1(t) % m*x2dd + k(2*x2-x1) + c(2*x2d - x1d) = f2(t) +% This example is described in Section 4.1 (Figures 1,2) of the following +% article +% S. Jain, T. Breunung, G. Haller, Fast Computation of Steady-State +% Response for Nonlinear Vibrations of High-Degree-of-Freedom Systems, +% Nonlinear Dyn (2019) 97: 313. https://doi.org/10.1007/s11071-019-04971-1 + + clear all clc; close all diff --git a/examples/ModifiedShawPierreQuasiPeriodic.m b/examples/ModifiedShawPierreQuasiPeriodic.m index a2346dd..02703d5 100644 --- a/examples/ModifiedShawPierreQuasiPeriodic.m +++ b/examples/ModifiedShawPierreQuasiPeriodic.m @@ -2,6 +2,13 @@ % m*x1dd + k(x1-x2) + c(x1d-x2d) + g*x1^3 = f1(t) % m*x2dd + k(2*x2-x1) + c(2*x2d - x1d) = f2(t) +% This example is described in Section 4.1.2 (Figure 4) of the following +% article +% S. Jain, T. Breunung, G. Haller, Fast Computation of Steady-State +% Response for Nonlinear Vibrations of High-Degree-of-Freedom Systems, +% Nonlinear Dyn (2019) 97: 313. https://doi.org/10.1007/s11071-019-04971-1 + + clear all clc; close all diff --git a/examples/motiv_example.m b/examples/motiv_example.m new file mode 100644 index 0000000..acad724 --- /dev/null +++ b/examples/motiv_example.m @@ -0,0 +1,101 @@ +%% Mass-spring motivational example for mode selection + +% This example is described in Section 4.3 (Figure 2) of the following +% article +% G. Buza, S. Jain, G. Haller, Using Spectral Submanifolds for Optimal +% Mode Selection in Model Reduction, (2020) Preprint available on arXiv.org + +clear +clc +close all + +%% properties +M = eye(3); +D1 = 0.01; +D2 = 0.02; +D3 = 0.08; +ome1 = 2; +ome2 = 3; +ome3 = 5; +K = diag([ome1^2 ome2^2 ome3^2]); +C = 2*diag([ome1*D1 ome2*D2 ome3*D3]); + +% Nonlinearity +S = @(x)SP_nonlinearity3(x,ome1,ome2,ome3); +DS = @(x)SP_nonlinearity_derv3(x,ome1,ome2,ome3); + +% Loading +A = 0.02; % loading amplitude +f0 = [1;0;0]; % loading shape +f = @(t,T)A*f0*sin(2*pi*t/T); % loading function + +%% Mode selection + +SS = SSR(M,C,K,1); + +% SSM computation +R = cell(1,2); +R{1} = ome2^2/2; +R{2} = ome3^2/2; +W = SS.SSM2(R); +disp(['norm(W_2)=' num2str(norm(W{1})) ', norm(W_3)=' num2str(norm(W{2}))]) +[~, ind] = max([norm(W{1}) norm(W{2})]); + +% master mode sets considered +I_1 = [1 2]; % mode selection 1 +I_2 = [1 ind+1]; % mode selection 2 + +%% SSR Package +SSfull = SSR(M,C,K,[1 2 3]); % full system +SSred1 = SSR(M,C,K,I_1); % reduced to I_1 +SSred2 = SSR(M,C,K,I_2); % reduced to I_2 +SSfull.S = S; +SSfull.DS = DS; +SSfull.f = f; +SSred1 = copyfromfull(SSred1,SSfull); +SSred2 = copyfromfull(SSred2,SSfull); + +%% Sequential continuation +omega1 = min(SSfull.omega); +Omega_range = [0.7 1.3]*omega1; +T_range = 2*pi./Omega_range; +ncontsteps = 100; + +[Omega_arrayf,Nf] = SSfull.FRC(T_range,ncontsteps); +[Omega_array1,N1] = SSred1.FRC(T_range,ncontsteps); +[Omega_array2,N2] = SSred2.FRC(T_range,ncontsteps); + +Nlin = nan(length(Omega_arrayf),1); +for j = 1:length(Omega_arrayf) + SSfull.T = 2*pi/Omega_arrayf(j); + [x_lin, xd_lin] = SSfull.LinearResponse(); + LinSol = [x_lin; xd_lin]; + Nlin(j) = SSfull.dt*norm(LinSol,'fro'); +end + +figure(2); semilogy(Omega_arrayf,Nf,'-k','DisplayName', 'Full','linewidth',1); axis tight; grid on; hold on; +xlabel('$$\Omega$$ [rad/s]'); ylabel('$$||q||_2$$') +legend('show') +semilogy(Omega_arrayf,Nlin,'-b', 'DisplayName', 'Linear','linewidth',1); +semilogy(Omega_array1,N1,'--g', 'DisplayName', 'Reduced $$I_1$$','linewidth',2); +semilogy(Omega_array2,N2,'--r','DisplayName', 'Reduced $$I_2$$','linewidth',2); + + +function f = SP_nonlinearity3(x,o1,o2,o3) +f = [o1^2/2*(3*x(1)^2+x(2)^2+x(3)^2)+o2^2*x(1)*x(2)+o3^2*x(1)*x(3)+(o1^2+o2^2+o3^2)/2*x(1)*(x(1)^2+x(2)^2+x(3)^2);... + o2^2/2*(3*x(2)^2+x(1)^2+x(3)^2)+o1^2*x(1)*x(2)+o3^2*x(2)*x(3)+(o1^2+o2^2+o3^2)/2*x(2)*(x(1)^2+x(2)^2+x(3)^2);... + o3^2/2*(3*x(3)^2+x(2)^2+x(1)^2)+o2^2*x(3)*x(2)+o1^2*x(1)*x(3)+(o1^2+o2^2+o3^2)/2*x(3)*(x(1)^2+x(2)^2+x(3)^2) ]; +end + +function Df = SP_nonlinearity_derv3(x,o1,o2,o3) +Df = zeros(3); +Df(1,1) = 3*o1^2*x(1)+o2^2*x(2)+o3^2*x(3) + (o1^2+o2^2+o3^2)/2*(3*x(1)^2 + x(2)^2+x(3)^2); +Df(2,2) = 3*o2^2*x(2)+o1^2*x(1)+o3^2*x(3) + (o1^2+o2^2+o3^2)/2*(3*x(2)^2 + x(1)^2+x(3)^2); +Df(3,3) = (3*o3^2*x(3)+o2^2*x(2)+o1^2*x(1) + (o1^2+o2^2+o3^2)/2*(3*x(3)^2 + x(2)^2+x(1)^2)); +Df(1,2) = o1^2*x(2)+o2^2*x(1)+(o1^2+o2^2+o3^2)*x(1)*x(2); +Df(1,3) = o1^2*x(3)+o3^2*x(1)+(o1^2+o2^2+o3^2)*x(1)*x(3); +Df(2,3) = o2^2*x(3)+o3^2*x(2)+(o1^2+o2^2+o3^2)*x(2)*x(3); +Df(3,2) = Df(2,3); +Df(2,1) = Df(1,2); +Df(3,1) = Df(1,3); +end \ No newline at end of file diff --git a/examples/ndof_oscillator_chain/OscillatorChain.m b/examples/ndof_oscillator_chain/OscillatorChain.m index d1df4fb..d6007fc 100644 --- a/examples/ndof_oscillator_chain/OscillatorChain.m +++ b/examples/ndof_oscillator_chain/OscillatorChain.m @@ -1,4 +1,10 @@ %% Oscillator chain with coupled nonlinearities +% This example is described in Section 4.2 (Figure 6) of the following +% article +% S. Jain, T. Breunung, G. Haller, Fast Computation of Steady-State +% Response for Nonlinear Vibrations of High-Degree-of-Freedom Systems, +% Nonlinear Dyn (2019) 97: 313. https://doi.org/10.1007/s11071-019-04971-1 + clear all clc diff --git a/examples/nonsmooth/NonSmoothOscillator.m b/examples/nonsmooth/NonSmoothOscillator.m index c81b7c3..dd6d739 100644 --- a/examples/nonsmooth/NonSmoothOscillator.m +++ b/examples/nonsmooth/NonSmoothOscillator.m @@ -1,4 +1,11 @@ %% 2-DOF oscillator with non smooth nonlinearity +% This example is described in Section 4.1.3 (Figure 5) of the following +% article +% S. Jain, T. Breunung, G. Haller, Fast Computation of Steady-State +% Response for Nonlinear Vibrations of High-Degree-of-Freedom Systems, +% Nonlinear Dyn (2019) 97: 313. https://doi.org/10.1007/s11071-019-04971-1 + + clear all clc; close all diff --git a/examples/vonKarman/example_curvedbeam.m b/examples/vonKarman/example_curvedbeam.m new file mode 100644 index 0000000..ca31e7a --- /dev/null +++ b/examples/vonKarman/example_curvedbeam.m @@ -0,0 +1,99 @@ +%% Mode Selection on Von Karman curved Beam Example +% This example is described in Section 6.2 (Figure 8) of the following +% article. +% G. Buza, S. Jain, G. Haller, Using Spectral Submanifolds for Optimal +% Mode Selection in Model Reduction, (2020) Preprint available on arXiv.org + +% The finite element model is taken from the following article. +% Jain, S., Tiso, P., & Haller, G. (2018). Exact nonlinear model reduction +% for a von Kármán beam: slow-fast decomposition and spectral submanifolds. +% Journal of Sound and Vibration, 423, 195–211. +% https://doi.org/10.1016/J.JSV.2018.01.049 + +clear +clc +close all + +%% parameters +epsilon = 7e-3; % beam thickness-to-length ratio +midp_height = 0.005; % height of midpoint relative to the ends (measure of curvature of the beam) +[Geo] = Geometry(epsilon,midp_height); + +nElements = 10; +BC = 'B'; % boundary condition: simply supported + +%% build finite element model +[Misc,model] = Beam_Model(Geo,nElements,BC); + +n = length(model.freeDOFs); % total number of degrees of freedom + +%% Mass Matrix +M = model.M(model.freeDOFs,model.freeDOFs); + +%% Stiffness Matrix +K = model.K(model.freeDOFs,model.freeDOFs); + +%% Damping Matrix +C = model.C(model.freeDOFs,model.freeDOFs); + +%% Nonlinearity +S =@(x) NonlinearityVK(model,x); +DS =@(x) NonlinearityJacobianVK(model,x); + +A = 80; % amplitude for which the FRC is computed +f = @(t,T) A*model.loads(model.freeDOFs) * sin(2*pi*t/T); % loading function + +%% Steady-State tool +SSfull = SSR(M,C,K,1:n); % full system + +%% Mode selection + +% Input parameters for mode selection +param.curvtolerance = 0.05; +param.nmodes = 10; +param.type = 1; % set to 0 if param.nmodes is the max number of modes, +% set to 1 if param.nmodes is the desired number of modes +param.example = 'simple'; % always use simple if enough memory is available to store X below +X = S11(model); % this conains the quadratic nonlinearities for all dofs (X_i). +% X_i is symmetric, and can be found from the form M\ddot{q} + C\dot{q} + +% Kq +\Gamma(q) = F(t) (i denotes the specific row of this equation) + +% master mode sets considered +modesel1 = 1:10; +modesel2 = modeselect(SSfull,param,X); +SSred1 = SSR(M,C,K,modesel1); % reduced to I_1 +SSred2 = SSR(M,C,K,modesel2); % reduced to I_2 + +SSfull.S = S; +SSfull.DS = DS; +SSfull.f = f; +SSred1 = copyfromfull(SSred1,SSfull); +SSred2 = copyfromfull(SSred2,SSfull); + +%% Sequential continuation +omega1 = SSfull.omega(1); +Omega_range = [0.7 1.3]*omega1; +T_range = 2*pi./Omega_range; +ncontsteps = 100; + +[Omega_arrayf,Nf] = SSfull.FRC(T_range,ncontsteps); +[Omega_array1,N1] = SSred1.FRC(T_range,ncontsteps); +[Omega_array2,N2] = SSred2.FRC(T_range,ncontsteps); + +Nlin = nan(length(Omega_arrayf),1); +for j = 1:length(Omega_arrayf) + SSfull.T = 2*pi/Omega_arrayf(j); + [x_lin, xd_lin] = SSfull.LinearResponse(); + LinSol = [x_lin; xd_lin]; + Nlin(j) = SSfull.dt*norm(LinSol,'fro'); +end + +figure(2); semilogy(Omega_arrayf,Nf,'-k','DisplayName', 'Full','linewidth',1); axis tight; grid on; hold on; +xlabel('$$\Omega$$ [rad/s]'); ylabel('$$||q||_2$$') +legend('show') +semilogy(Omega_arrayf,Nlin,'-b', 'DisplayName', 'Linear','linewidth',1); +semilogy(Omega_array1,N1,'--g', 'DisplayName', 'Reduced $$I_1$$','linewidth',2); +semilogy(Omega_array2,N2,'--r','DisplayName', 'Reduced $$I_2$$','linewidth',2); + + + diff --git a/examples/vonKarman/example_straightbeam.m b/examples/vonKarman/example_straightbeam.m new file mode 100644 index 0000000..25c298f --- /dev/null +++ b/examples/vonKarman/example_straightbeam.m @@ -0,0 +1,136 @@ +%% Mode selection on Von Karman Straight Beam example +% This example is described in Section 6.1 (Figure 5) of the following +% article. +% G. Buza, S. Jain, G. Haller, Using Spectral Submanifolds for Optimal +% Mode Selection in Model Reduction, (2020) Preprint available on arXiv.org + +% The finite element model is taken from the following article. +% Jain, S., Tiso, P., & Haller, G. (2018). Exact nonlinear model reduction +% for a von Kármán beam: slow-fast decomposition and spectral submanifolds. +% Journal of Sound and Vibration, 423, 195–211. +% https://doi.org/10.1016/J.JSV.2018.01.049 + +clear +clc +close all +%% parameters +epsilon = 1e-3; % beam thickness-to-length ratio +midp_height = 0; % height of midpoint relative to the ends (measure of curvature of the beam) +[Geo] = Geometry(epsilon,midp_height); + +nElements = 10; +BC = 'B'; % boundary condition: simply supported + +%% build finite element model +[Misc,model] = Beam_Model(Geo,nElements,BC); + +n = length(model.freeDOFs); % total number of degrees of freedom +plot_dof = 9; % the degree of freedom at which the response is plotted +%% Mass Matrix +M = model.M(model.freeDOFs,model.freeDOFs); + +%% Stiffness Matrix +K = model.K(model.freeDOFs,model.freeDOFs); + +%% Damping Matrix +C = model.C(model.freeDOFs,model.freeDOFs); + +%% Steady-State tool +SS = SSR(M,C,K,1:n); % full system + +%% Mode selection + +% Input parameters for mode selection +param.curvtolerance = 0.05; +param.nmodes = 10; +param.type = 0; % set to 0 if param.nmodes is the max number of modes, + % set to 1 if param.nmodes is the desired number of modes +param.example = 'simple'; % always use simple if enough memory is available to store X below +param.initialmodes = [1 2 3 4 5]; % optional input modes +X = S11(model); % this conains the quadratic nonlinearities for all dofs (X_i). + % X_i is symmetric, and can be found from the form M\ddot{q} + C\dot{q} + + % Kq +\Gamma(q) = F(t) (i denotes the specific row of this equation) + +% master mode sets considered +I_1 = 1:10; +I_2 = modeselect(SS,param,X); +SS1 = SSR(M,C,K,I_1); % reduced to I_1 +SS2 = SSR(M,C,K,I_2); % reduced to I_2 + +%% Incrementing amplitude +T = 2*pi/26; +A = [10:5:35 40:4:172 174 176:1:200 200:1:230]* 10^(-2); % amplitude vector +SS.S = @(x) NonlinearityVK(model,x); +SS.DS = @(x) NonlinearityJacobianVK(model,x); + +tol = 1e-6; % Newton Raphson error tolerance +maxiter = 60; +N2 = nan(length(A),1); +N1 = nan(length(A),1); +NF = nan(length(A),1); +NL = nan(length(A),1); + +for j = 1:length(A) + Fext = @(t,T) A(j)*model.loads(model.freeDOFs) * sin(2*pi*t/T); + SS.f = Fext; + SS.T = T; + SS.n_steps = 100; + SS1 = copyfromfull(SS1,SS); + SS2 = copyfromfull(SS2,SS); + + if j == 1 + [x2,xd2] = SS1.NewtonRaphson('tol',tol,'maxiter',maxiter); % find out response using Picard iteration + else + [x2,xd2] = SS1.NewtonRaphson('init',x2,'tol',tol,'maxiter',maxiter); + end + N1(j) = SS1.dt*norm([x2; xd2],'fro'); + + if j == 1 + [x3,xd3] = SS2.NewtonRaphson('tol',tol,'maxiter',maxiter); % find out response using Picard iteration + else + [x3,xd3] = SS2.NewtonRaphson('init',x3,'tol',tol,'maxiter',maxiter); + end + N2(j) = SS2.dt*norm([x3; xd3],'fro'); + + + if j == 1 + [xf,xdf] = SS.NewtonRaphson('tol',tol,'maxiter',maxiter); + else + [xf,xdf] = SS.NewtonRaphson('init',xf,'tol',tol,'maxiter',maxiter); + end + + NF(j) = SS.dt*norm([xf; xdf],'fro'); + + [x_lin, xd_lin] = SS.LinearResponse(); + NL(j) = SS.dt*norm([x_lin; xd_lin],'fro'); + + disp(vpa(A(j))) +end + +figure(3); + +%% final amplitude time history +plot(SS.t, xf(plot_dof,:), 'k', 'DisplayName', 'Full','linewidth',1); +axis tight; +legend('show') +ylim([-3.1*10^(-3) 3.1*10^(-3)]) +xlabel('time (seconds)'); +ylabel('displacement $$w_{4}$$ [m]') +hold on; +plot(SS.t, x_lin(plot_dof,:), 'b', 'DisplayName', 'Linear','linewidth',1); +plot(SS1.t, x2(plot_dof,:),'--g', 'DisplayName', 'Reduced $$I_1$$','linewidth',2); +plot(SS2.t, x3(plot_dof,:),'--r', 'DisplayName', 'Reduced $$I_2$$','linewidth',2); + +%% forcing amp - response amp curve +figure +plot(A,NF,'-k','DisplayName', 'Full','linewidth',1); axis tight; hold on; +xlabel('$$F$$ [N]'); ylabel('$$||q||_2$$') +lgd = legend('show'); +lgd.Location = 'southeast'; +plot(A,NL,'-b', 'DisplayName', 'Linear','linewidth',1); +plot(A,N1,'--g', 'DisplayName', 'Reduced $$I_1$$','linewidth',2); +plot(A,N2,'--r','DisplayName', 'Reduced $$I_2$$','linewidth',2); +ylim([0 0.032]) + + + diff --git a/examples/vonKarman/model/AssembleLinearStiffness.m b/examples/vonKarman/model/AssembleLinearStiffness.m new file mode 100644 index 0000000..d27a5b0 --- /dev/null +++ b/examples/vonKarman/model/AssembleLinearStiffness.m @@ -0,0 +1,33 @@ +function [K] = AssembleLinearStiffness(model) +% The finite element model is taken from the following article: +% Jain, S., Tiso, P., & Haller, G. (2018). Exact nonlinear model reduction +% for a von Kármán beam: slow-fast decomposition and spectral submanifolds. +% Journal of Sound and Vibration, 423, 195–211. +% https://doi.org/10.1016/J.JSV.2018.01.049 + +ndof = model.nDOF; +nelem = model.nElements; +I = zeros(nelem*6*6,1); +J = zeros(nelem*6*6,1); +KK = zeros(nelem*6*6,1); + +el = 0; +% F = zeros(ndof,1); +for i = 1:model.nElements + index = 3*(i-1)+1:3*(i+1); + + Te = compute_rotation_matrix(model.nodes(i:i+1,:)); + + [Kel,~] = ElementStiffness(model,zeros(6,1)); + Kel = Te.' * Kel * Te; + +% F(index) = F(index)+ Fel; + for col = 1:6 + I(el+6*(col-1)+1:el+6*col,1) = ones(6,1)*index(col); + J(el+6*(col-1)+1:el+6*col,1) = index; + KK(el+6*(col-1)+1:el+6*col,1) = Kel(col,:); + end + el = el + 36; + +end +K = sparse(I,J,KK,ndof,ndof); diff --git a/examples/vonKarman/model/AssembleMass.m b/examples/vonKarman/model/AssembleMass.m new file mode 100644 index 0000000..607fbd7 --- /dev/null +++ b/examples/vonKarman/model/AssembleMass.m @@ -0,0 +1,30 @@ +function M = AssembleMass(model) +% The finite element model is taken from the following article: +% Jain, S., Tiso, P., & Haller, G. (2018). Exact nonlinear model reduction +% for a von Kármán beam: slow-fast decomposition and spectral submanifolds. +% Journal of Sound and Vibration, 423, 195–211. +% https://doi.org/10.1016/J.JSV.2018.01.049 + +ndof = model.nDOF; +nelem = model.nElements; +I = zeros(nelem*6*6,1); +J = zeros(nelem*6*6,1); +KK = zeros(nelem*6*6,1); +el = 0; +for i = 1:model.nElements + index = 3*(i-1)+1:3*(i+1); + + Te = compute_rotation_matrix(model.nodes(i:i+1,:)); + + Mel = ElementMass(model); + Mel = Te.' * Mel * Te; + + for col = 1:6 + I(el+6*(col-1)+1:el+6*col,1) = ones(6,1)*index(col); + J(el+6*(col-1)+1:el+6*col,1) = index; + KK(el+6*(col-1)+1:el+6*col,1) = Mel(col,:); + end + el = el + 36; + +end +M = sparse(I,J,KK,ndof,ndof); diff --git a/examples/vonKarman/model/AssembleTangentStiffness.m b/examples/vonKarman/model/AssembleTangentStiffness.m new file mode 100644 index 0000000..dcd42e3 --- /dev/null +++ b/examples/vonKarman/model/AssembleTangentStiffness.m @@ -0,0 +1,35 @@ +function [K,F] = AssembleTangentStiffness(model,u) +% The finite element model is taken from the following article: +% Jain, S., Tiso, P., & Haller, G. (2018). Exact nonlinear model reduction +% for a von Kármán beam: slow-fast decomposition and spectral submanifolds. +% Journal of Sound and Vibration, 423, 195–211. +% https://doi.org/10.1016/J.JSV.2018.01.049 + +ndof = model.nDOF; +nelem = model.nElements; +I = zeros(nelem*6*6,1); +J = zeros(nelem*6*6,1); +KK = zeros(nelem*6*6,1); + +el = 0; +F = zeros(ndof,1); +for i = 1:model.nElements + index = 3*(i-1)+1:3*(i+1); + + Te = compute_rotation_matrix(model.nodes(i:i+1,:)); + + qe = Te * u(index); + [Kel,Fel] = ElementStiffness(model,qe); + + Kel = Te.' * Kel * Te; + F(index) = F(index)+ Te.' * Fel; + + for col = 1:6 + I(el+6*(col-1)+1:el+6*col,1) = ones(6,1)*index(col); + J(el+6*(col-1)+1:el+6*col,1) = index; + KK(el+6*(col-1)+1:el+6*col,1) = Kel(col,:); + end + el = el + 36; + +end +K = sparse(I,J,KK,ndof,ndof); diff --git a/examples/vonKarman/model/Beam_Model.m b/examples/vonKarman/model/Beam_Model.m new file mode 100644 index 0000000..3ca72f9 --- /dev/null +++ b/examples/vonKarman/model/Beam_Model.m @@ -0,0 +1,89 @@ +function [Misc,model] = Beam_Model(Geo,nElements,B) +% This finite element model is taken from the following article: +% Jain, S., Tiso, P., & Haller, G. (2018). Exact nonlinear model reduction +% for a von Kármán beam: slow-fast decomposition and spectral submanifolds. +% Journal of Sound and Vibration, 423, 195–211. +% https://doi.org/10.1016/J.JSV.2018.01.049 + +Damping = 2; % Structural damping : 0 - Nil, 1 - Rayleigh, 2 - Linear viscoelastic + +% Assumed a curved beam with N equidistant elements along the x axis +model.nElements = nElements; +model.nNodes = nElements +1; +model.nDOF = model.nNodes*3; +model.dx = Geo.l/model.nElements; +model.b = Geo.b; + +model.nodes = [(0:model.dx:Geo.l).' zeros(model.nNodes,1)]; +if Geo.a~=0 + R = 1/2*(Geo.a^2+(Geo.l/2)^2)/Geo.a; + theta0 = asin(Geo.l/2/R); + + for i = 1:size(model.nodes,1) + th = asin((model.nodes(i,1)-Geo.l/2)/R); + model.nodes(i,2) = R * cos(th) - R * cos(theta0); + end +end + +% Material Properties (Aluminium) +model.E = 70e9; % Youngs Modulus (Pa) +model.A = Geo.A; % Cross section area of Beam +model.I = Geo.I; % Moment of inertia of beam +model.rho = 2700; % Density +model.kappa = 1e8; % Material damping for Kelvin Voigt model + + +% Loading +model.pw = 1; % loading factor along the transverse direction +model.pu = 0; % loading factor along the axial direction + +model.loads = sparse(model.nDOF,1); +l = ones(model.nNodes,1); +l(2:end-1) = 2; +model.uDOF = 1:3:model.nDOF; % axial degrees of freedom +model.wDOF = setdiff(1:model.nDOF,model.uDOF); % transverse DOFs +model.loads(2:3:end) = model.pw*(model.dx/2)*l; % +model.loads(model.uDOF) = model.pu*(model.dx/2)*l; + +% Boundary conditions +model.fixDOFs = BoundaryCoundition(model,B); +model.freeDOFs = setdiff(1:model.nDOF, model.fixDOFs); +model.uDOF = setdiff(model.uDOF,model.fixDOFs); +model.wDOF = setdiff(model.wDOF,model.fixDOFs); + +% Stiffness +model.K = AssembleLinearStiffness(model); +model.M = AssembleMass(model); + +% +model.modes = 1:5; +model.m = length(model.modes); +% Vibration problem +Mred = model.M(model.freeDOFs,model.freeDOFs); +model.Ks = model.K; +Kred = model.K(model.freeDOFs,model.freeDOFs); +[model.V,D] = eigs(Kred,Mred,max(model.modes),'SM'); % Return M smallest eigenvalues and vectors +[omega2,V] = sortVMsVK(model.V,Mred,D); +disp('First 5 undamped Natural Frequencies [Hz]') +disp(sqrt(omega2)/(2*pi)) +model.omega2 = omega2(model.modes); +model.VMs = zeros(model.nDOF,model.m); +model.VMs(model.freeDOFs,:) = V(:,model.modes); +Misc.omega = sqrt(model.omega2); + + +figure +PlotDeformation(model,model.VMs(:,1),0.1); + +switch Damping + case 1 + % Rayleigh Damping + W = sqrt(model.omega2(1:2)); + a = [W(1) 1/W(1);W(2) 1/W(2)]\[0.02;0.02]; + model.C = (a(2)*model.M + a(1)*model.K); + + case 2 + % ViscoElastic Damping with linear strain rate + model.C = (model.kappa/model.E)*model.Ks; +end + diff --git a/examples/vonKarman/model/BoundaryCoundition.m b/examples/vonKarman/model/BoundaryCoundition.m new file mode 100644 index 0000000..4fa2e49 --- /dev/null +++ b/examples/vonKarman/model/BoundaryCoundition.m @@ -0,0 +1,11 @@ +function [fixDOFs] = BoundaryCoundition(model,S) + +switch S + case 'C' % Cantilever + fixDOFs = 1:3; + case 'B' % Simply supported + fixDOFs = [1 2 model.nDOF-2 model.nDOF-1]; + case 'DC' % Doubly clamped + fixDOFs = [1:3 (model.nDOF-2):model.nDOF]; +end + \ No newline at end of file diff --git a/examples/vonKarman/model/ElementMass.m b/examples/vonKarman/model/ElementMass.m new file mode 100644 index 0000000..ce6a432 --- /dev/null +++ b/examples/vonKarman/model/ElementMass.m @@ -0,0 +1,22 @@ +function M = ElementMass(prop) +% The finite element model is taken from the following article: +% Jain, S., Tiso, P., & Haller, G. (2018). Exact nonlinear model reduction +% for a von Kármán beam: slow-fast decomposition and spectral submanifolds. +% Journal of Sound and Vibration, 423, 195–211. +% https://doi.org/10.1016/J.JSV.2018.01.049 + +rho = prop.rho; +A = prop.A; +l = prop.dx; + + +M = sparse(6,6); + +M([1 4],[1 4]) = ((rho*A*l)/6)*[2 1; + 1 2]; + +M([2 3 5 6],[2 3 5 6]) = ((rho*A*l)/420)*[156 22*l 54 -13*l; + 22*l 4*l^2 13*l -3*l^2; + 54 13*l 156 -22*l; + -13*l -3*l^2 -22*l 4*l^2]; + diff --git a/examples/vonKarman/model/ElementStiffness.m b/examples/vonKarman/model/ElementStiffness.m new file mode 100644 index 0000000..b6ede00 --- /dev/null +++ b/examples/vonKarman/model/ElementStiffness.m @@ -0,0 +1,34 @@ +function [K,F] = ElementStiffness(model,q) +% The finite element model is taken from the following article: +% Jain, S., Tiso, P., & Haller, G. (2018). Exact nonlinear model reduction +% for a von Kármán beam: slow-fast decomposition and spectral submanifolds. +% Journal of Sound and Vibration, 423, 195–211. +% https://doi.org/10.1016/J.JSV.2018.01.049 + +E = model.E; +A = model.A; +dx = model.dx; +I = model.I; + +u1 = q(1); +w1 = q(2); +t1 = q(3); +u2 = q(4); +w2 = q(5); +t2 = q(6); + + +K = [ (A*E)/dx, -((A*E*(36*w1 - 36*w2))/30 + (A*E*dx*(3*t1 + 3*t2))/30)/dx^2, - (A*E*(4*t1 - t2))/30 - (A*E*(3*w1 - 3*w2))/(30*dx), -(A*E)/dx, ((A*E*(36*w1 - 36*w2))/30 + (A*E*dx*(3*t1 + 3*t2))/30)/dx^2, (A*E*(t1 - 4*t2))/30 - (A*E*(3*w1 - 3*w2))/(30*dx); + (6*A*E*w2)/(5*dx^2) - (A*E*t2)/(10*dx) - (6*A*E*w1)/(5*dx^2) - (A*E*t1)/(10*dx), (12*E*I)/dx^3 - (6*A*E*u1)/(5*dx^2) + (6*A*E*u2)/(5*dx^2) + (9*A*E*t1^2)/(70*dx) + (9*A*E*t2^2)/(70*dx) + (108*A*E*w1^2)/(35*dx^3) + (108*A*E*w2^2)/(35*dx^3) + (27*A*E*t1*w1)/(35*dx^2) - (27*A*E*t1*w2)/(35*dx^2) + (27*A*E*t2*w1)/(35*dx^2) - (27*A*E*t2*w2)/(35*dx^2) - (216*A*E*w1*w2)/(35*dx^3), (6*E*I)/dx^2 - (3*A*E*t1^2)/280 + (3*A*E*t2^2)/280 - (A*E*u1)/(10*dx) + (A*E*u2)/(10*dx) + (27*A*E*w1^2)/(70*dx^2) + (27*A*E*w2^2)/(70*dx^2) + (3*A*E*t1*t2)/140 + (9*A*E*t1*w1)/(35*dx) - (9*A*E*t1*w2)/(35*dx) - (27*A*E*w1*w2)/(35*dx^2), (A*E*t1)/(10*dx) + (A*E*t2)/(10*dx) + (6*A*E*w1)/(5*dx^2) - (6*A*E*w2)/(5*dx^2), (6*A*E*u1)/(5*dx^2) - (12*E*I)/dx^3 - (6*A*E*u2)/(5*dx^2) - (9*A*E*t1^2)/(70*dx) - (9*A*E*t2^2)/(70*dx) - (108*A*E*w1^2)/(35*dx^3) - (108*A*E*w2^2)/(35*dx^3) - (27*A*E*t1*w1)/(35*dx^2) + (27*A*E*t1*w2)/(35*dx^2) - (27*A*E*t2*w1)/(35*dx^2) + (27*A*E*t2*w2)/(35*dx^2) + (216*A*E*w1*w2)/(35*dx^3), (6*E*I)/dx^2 + (3*A*E*t1^2)/280 - (3*A*E*t2^2)/280 - (A*E*u1)/(10*dx) + (A*E*u2)/(10*dx) + (27*A*E*w1^2)/(70*dx^2) + (27*A*E*w2^2)/(70*dx^2) + (3*A*E*t1*t2)/140 + (9*A*E*t2*w1)/(35*dx) - (9*A*E*t2*w2)/(35*dx) - (27*A*E*w1*w2)/(35*dx^2); + (A*E*t2)/30 - (2*A*E*t1)/15 - (A*E*w1)/(10*dx) + (A*E*w2)/(10*dx), (6*E*I)/dx^2 - (3*A*E*t1^2)/280 + (3*A*E*t2^2)/280 - (A*E*u1)/(10*dx) + (A*E*u2)/(10*dx) + (27*A*E*w1^2)/(70*dx^2) + (27*A*E*w2^2)/(70*dx^2) + (3*A*E*t1*t2)/140 + (9*A*E*t1*w1)/(35*dx) - (9*A*E*t1*w2)/(35*dx) - (27*A*E*w1*w2)/(35*dx^2), (3*A*E*dx*t1^2)/35 - (3*A*E*dx*t1*t2)/140 - (3*A*E*t1*w1)/140 + (3*A*E*t1*w2)/140 + (A*E*dx*t2^2)/140 + (3*A*E*t2*w1)/140 - (3*A*E*t2*w2)/140 + (9*A*E*w1^2)/(70*dx) - (9*A*E*w1*w2)/(35*dx) + (9*A*E*w2^2)/(70*dx) - (2*A*E*u1)/15 + (2*A*E*u2)/15 + (4*E*I)/dx, (2*A*E*t1)/15 - (A*E*t2)/30 + (A*E*w1)/(10*dx) - (A*E*w2)/(10*dx), (3*A*E*t1^2)/280 - (6*E*I)/dx^2 - (3*A*E*t2^2)/280 + (A*E*u1)/(10*dx) - (A*E*u2)/(10*dx) - (27*A*E*w1^2)/(70*dx^2) - (27*A*E*w2^2)/(70*dx^2) - (3*A*E*t1*t2)/140 - (9*A*E*t1*w1)/(35*dx) + (9*A*E*t1*w2)/(35*dx) + (27*A*E*w1*w2)/(35*dx^2), (A*E*u1)/30 - (A*E*u2)/30 + (2*E*I)/dx - (3*A*E*dx*t1^2)/280 - (3*A*E*dx*t2^2)/280 + (3*A*E*t1*w1)/140 - (3*A*E*t1*w2)/140 + (3*A*E*t2*w1)/140 - (3*A*E*t2*w2)/140 + (A*E*dx*t1*t2)/70; + -(A*E)/dx, ((A*E*(36*w1 - 36*w2))/30 + (A*E*dx*(3*t1 + 3*t2))/30)/dx^2, (A*E*(4*t1 - t2))/30 + (A*E*(3*w1 - 3*w2))/(30*dx), (A*E)/dx, -((A*E*(36*w1 - 36*w2))/30 + (A*E*dx*(3*t1 + 3*t2))/30)/dx^2, (A*E*(3*w1 - 3*w2))/(30*dx) - (A*E*(t1 - 4*t2))/30; + (A*E*t1)/(10*dx) + (A*E*t2)/(10*dx) + (6*A*E*w1)/(5*dx^2) - (6*A*E*w2)/(5*dx^2), (6*A*E*u1)/(5*dx^2) - (12*E*I)/dx^3 - (6*A*E*u2)/(5*dx^2) - (9*A*E*t1^2)/(70*dx) - (9*A*E*t2^2)/(70*dx) - (108*A*E*w1^2)/(35*dx^3) - (108*A*E*w2^2)/(35*dx^3) - (27*A*E*t1*w1)/(35*dx^2) + (27*A*E*t1*w2)/(35*dx^2) - (27*A*E*t2*w1)/(35*dx^2) + (27*A*E*t2*w2)/(35*dx^2) + (216*A*E*w1*w2)/(35*dx^3), (3*A*E*t1^2)/280 - (6*E*I)/dx^2 - (3*A*E*t2^2)/280 + (A*E*u1)/(10*dx) - (A*E*u2)/(10*dx) - (27*A*E*w1^2)/(70*dx^2) - (27*A*E*w2^2)/(70*dx^2) - (3*A*E*t1*t2)/140 - (9*A*E*t1*w1)/(35*dx) + (9*A*E*t1*w2)/(35*dx) + (27*A*E*w1*w2)/(35*dx^2), (6*A*E*w2)/(5*dx^2) - (A*E*t2)/(10*dx) - (6*A*E*w1)/(5*dx^2) - (A*E*t1)/(10*dx), (12*E*I)/dx^3 - (6*A*E*u1)/(5*dx^2) + (6*A*E*u2)/(5*dx^2) + (9*A*E*t1^2)/(70*dx) + (9*A*E*t2^2)/(70*dx) + (108*A*E*w1^2)/(35*dx^3) + (108*A*E*w2^2)/(35*dx^3) + (27*A*E*t1*w1)/(35*dx^2) - (27*A*E*t1*w2)/(35*dx^2) + (27*A*E*t2*w1)/(35*dx^2) - (27*A*E*t2*w2)/(35*dx^2) - (216*A*E*w1*w2)/(35*dx^3), (3*A*E*t2^2)/280 - (3*A*E*t1^2)/280 - (6*E*I)/dx^2 + (A*E*u1)/(10*dx) - (A*E*u2)/(10*dx) - (27*A*E*w1^2)/(70*dx^2) - (27*A*E*w2^2)/(70*dx^2) - (3*A*E*t1*t2)/140 - (9*A*E*t2*w1)/(35*dx) + (9*A*E*t2*w2)/(35*dx) + (27*A*E*w1*w2)/(35*dx^2); + (A*E*t1)/30 - (2*A*E*t2)/15 - (A*E*w1)/(10*dx) + (A*E*w2)/(10*dx), (6*E*I)/dx^2 + (3*A*E*t1^2)/280 - (3*A*E*t2^2)/280 - (A*E*u1)/(10*dx) + (A*E*u2)/(10*dx) + (27*A*E*w1^2)/(70*dx^2) + (27*A*E*w2^2)/(70*dx^2) + (3*A*E*t1*t2)/140 + (9*A*E*t2*w1)/(35*dx) - (9*A*E*t2*w2)/(35*dx) - (27*A*E*w1*w2)/(35*dx^2), (A*E*u1)/30 - (A*E*u2)/30 + (2*E*I)/dx - (3*A*E*dx*t1^2)/280 - (3*A*E*dx*t2^2)/280 + (3*A*E*t1*w1)/140 - (3*A*E*t1*w2)/140 + (3*A*E*t2*w1)/140 - (3*A*E*t2*w2)/140 + (A*E*dx*t1*t2)/70, (2*A*E*t2)/15 - (A*E*t1)/30 + (A*E*w1)/(10*dx) - (A*E*w2)/(10*dx), (3*A*E*t2^2)/280 - (3*A*E*t1^2)/280 - (6*E*I)/dx^2 + (A*E*u1)/(10*dx) - (A*E*u2)/(10*dx) - (27*A*E*w1^2)/(70*dx^2) - (27*A*E*w2^2)/(70*dx^2) - (3*A*E*t1*t2)/140 - (9*A*E*t2*w1)/(35*dx) + (9*A*E*t2*w2)/(35*dx) + (27*A*E*w1*w2)/(35*dx^2), (A*E*dx*t1^2)/140 - (3*A*E*dx*t1*t2)/140 + (3*A*E*t1*w1)/140 - (3*A*E*t1*w2)/140 + (3*A*E*dx*t2^2)/35 - (3*A*E*t2*w1)/140 + (3*A*E*t2*w2)/140 + (9*A*E*w1^2)/(70*dx) - (9*A*E*w1*w2)/(35*dx) + (9*A*E*w2^2)/(70*dx) - (2*A*E*u1)/15 + (2*A*E*u2)/15 + (4*E*I)/dx]; + +F = [ - ((A*E*(18*w1^2 - 36*w1*w2 + 18*w2^2))/30 - (A*E*dx*(30*u1 - 30*u2 - 3*t1*w1 + 3*t1*w2 - 3*t2*w1 + 3*t2*w2))/30)/dx^2 - (A*E*(2*t1^2 - t1*t2 + 2*t2^2))/30; + (E*(3360*I*w1 - 3360*I*w2 + 288*A*w1^3 - 288*A*w2^3 - A*dx^3*t1^3 - A*dx^3*t2^3 + 1680*I*dx*t1 + 1680*I*dx*t2 + 864*A*w1*w2^2 - 864*A*w1^2*w2 - 28*A*dx^2*t1*u1 + 28*A*dx^2*t1*u2 - 28*A*dx^2*t2*u1 + 28*A*dx^2*t2*u2 + 108*A*dx*t1*w1^2 + 108*A*dx*t1*w2^2 + 108*A*dx*t2*w1^2 + 108*A*dx*t2*w2^2 + 3*A*dx^3*t1*t2^2 + 3*A*dx^3*t1^2*t2 + 36*A*dx^2*t1^2*w1 - 36*A*dx^2*t1^2*w2 + 36*A*dx^2*t2^2*w1 - 36*A*dx^2*t2^2*w2 - 336*A*dx*u1*w1 + 336*A*dx*u1*w2 + 336*A*dx*u2*w1 - 336*A*dx*u2*w2 - 216*A*dx*t1*w1*w2 - 216*A*dx*t2*w1*w2))/(280*dx^3); + (E*(5040*I*w1 - 5040*I*w2 + 108*A*w1^3 - 108*A*w2^3 + 24*A*dx^3*t1^3 - 3*A*dx^3*t2^3 + 3360*I*dx*t1 + 1680*I*dx*t2 + 324*A*w1*w2^2 - 324*A*w1^2*w2 - 112*A*dx^2*t1*u1 + 112*A*dx^2*t1*u2 + 28*A*dx^2*t2*u1 - 28*A*dx^2*t2*u2 + 108*A*dx*t1*w1^2 + 108*A*dx*t1*w2^2 + 6*A*dx^3*t1*t2^2 - 9*A*dx^3*t1^2*t2 - 9*A*dx^2*t1^2*w1 + 9*A*dx^2*t1^2*w2 + 9*A*dx^2*t2^2*w1 - 9*A*dx^2*t2^2*w2 - 84*A*dx*u1*w1 + 84*A*dx*u1*w2 + 84*A*dx*u2*w1 - 84*A*dx*u2*w2 - 216*A*dx*t1*w1*w2 + 18*A*dx^2*t1*t2*w1 - 18*A*dx^2*t1*t2*w2))/(840*dx^2); + ((A*E*(18*w1^2 - 36*w1*w2 + 18*w2^2))/30 - (A*E*dx*(30*u1 - 30*u2 - 3*t1*w1 + 3*t1*w2 - 3*t2*w1 + 3*t2*w2))/30)/dx^2 + (A*E*(2*t1^2 - t1*t2 + 2*t2^2))/30; + -(E*(3360*I*w1 - 3360*I*w2 + 288*A*w1^3 - 288*A*w2^3 - A*dx^3*t1^3 - A*dx^3*t2^3 + 1680*I*dx*t1 + 1680*I*dx*t2 + 864*A*w1*w2^2 - 864*A*w1^2*w2 - 28*A*dx^2*t1*u1 + 28*A*dx^2*t1*u2 - 28*A*dx^2*t2*u1 + 28*A*dx^2*t2*u2 + 108*A*dx*t1*w1^2 + 108*A*dx*t1*w2^2 + 108*A*dx*t2*w1^2 + 108*A*dx*t2*w2^2 + 3*A*dx^3*t1*t2^2 + 3*A*dx^3*t1^2*t2 + 36*A*dx^2*t1^2*w1 - 36*A*dx^2*t1^2*w2 + 36*A*dx^2*t2^2*w1 - 36*A*dx^2*t2^2*w2 - 336*A*dx*u1*w1 + 336*A*dx*u1*w2 + 336*A*dx*u2*w1 - 336*A*dx*u2*w2 - 216*A*dx*t1*w1*w2 - 216*A*dx*t2*w1*w2))/(280*dx^3); + (E*(5040*I*w1 - 5040*I*w2 + 108*A*w1^3 - 108*A*w2^3 - 3*A*dx^3*t1^3 + 24*A*dx^3*t2^3 + 1680*I*dx*t1 + 3360*I*dx*t2 + 324*A*w1*w2^2 - 324*A*w1^2*w2 + 28*A*dx^2*t1*u1 - 28*A*dx^2*t1*u2 - 112*A*dx^2*t2*u1 + 112*A*dx^2*t2*u2 + 108*A*dx*t2*w1^2 + 108*A*dx*t2*w2^2 - 9*A*dx^3*t1*t2^2 + 6*A*dx^3*t1^2*t2 + 9*A*dx^2*t1^2*w1 - 9*A*dx^2*t1^2*w2 - 9*A*dx^2*t2^2*w1 + 9*A*dx^2*t2^2*w2 - 84*A*dx*u1*w1 + 84*A*dx*u1*w2 + 84*A*dx*u2*w1 - 84*A*dx*u2*w2 - 216*A*dx*t2*w1*w2 + 18*A*dx^2*t1*t2*w1 - 18*A*dx^2*t1*t2*w2))/(840*dx^2)]; + diff --git a/examples/vonKarman/model/Geometry.m b/examples/vonKarman/model/Geometry.m new file mode 100644 index 0000000..dadfb23 --- /dev/null +++ b/examples/vonKarman/model/Geometry.m @@ -0,0 +1,9 @@ +function [Geo] = Geometry(eps,a) +%% Geometry for a Von Karman Beam element aligned along the x axis (all units SI) +Geo.eps = eps; +Geo.l = 1 ; % length +Geo.h = eps*Geo.l; % height +Geo.b = Geo.l/10; % width +Geo.A = Geo.b*Geo.h; % cross section area +Geo.I = Geo.b*Geo.h^3/12; % second moment of cross section +Geo.a = a; % height of midpoint relative to the ends (measure of curvature of the beam) diff --git a/examples/vonKarman/model/NonlinearityJacobianVK.m b/examples/vonKarman/model/NonlinearityJacobianVK.m new file mode 100644 index 0000000..915cea7 --- /dev/null +++ b/examples/vonKarman/model/NonlinearityJacobianVK.m @@ -0,0 +1,16 @@ +function [DS] = NonlinearityJacobianVK(model,x) +% The finite element model is taken from the following article: +% Jain, S., Tiso, P., & Haller, G. (2018). Exact nonlinear model reduction +% for a von Kármán beam: slow-fast decomposition and spectral submanifolds. +% Journal of Sound and Vibration, 423, 195–211. +% https://doi.org/10.1016/J.JSV.2018.01.049 + +q = zeros(model.nDOF,1); +q(model.freeDOFs) = x; + +[Kt,~] = AssembleTangentStiffness(model,q); + +DS = Kt - model.K; +DS = DS(model.freeDOFs,model.freeDOFs); +end + diff --git a/examples/vonKarman/model/NonlinearityVK.m b/examples/vonKarman/model/NonlinearityVK.m new file mode 100644 index 0000000..166ec62 --- /dev/null +++ b/examples/vonKarman/model/NonlinearityVK.m @@ -0,0 +1,14 @@ +function [S] = NonlinearityVK(model,x) +% The finite element model is taken from the following article: +% Jain, S., Tiso, P., & Haller, G. (2018). Exact nonlinear model reduction +% for a von Kármán beam: slow-fast decomposition and spectral submanifolds. +% Journal of Sound and Vibration, 423, 195–211. +% https://doi.org/10.1016/J.JSV.2018.01.049 + +q = zeros(model.nDOF,1); +q(model.freeDOFs) = x; +[~,F] = AssembleTangentStiffness(model,q); +S = F - model.K * q; +S = S(model.freeDOFs); +end + diff --git a/examples/vonKarman/model/PlotDeformation.m b/examples/vonKarman/model/PlotDeformation.m new file mode 100644 index 0000000..ef9fbd2 --- /dev/null +++ b/examples/vonKarman/model/PlotDeformation.m @@ -0,0 +1,25 @@ +function PlotDeformation(model,displacement,scale) + +l = model.nElements*model.dx; +u = displacement(1:3:end,:); w = displacement(2:3:end,:); +dis = sqrt(u.^2 + w.^2); +maxdispl = max(dis); + +sf = scale*l./maxdispl; +sf = diag(sf); +if strcmp(scale, 'TRUE') + sf = 1; +end +% sf = scale; + +x = 0:model.dx:l; +% y = zeros(size(x)); +xd = x' + u*sf; +yd = w*sf; + +% plot(x,y,'-k','LineWidth',2) +% hold on +plot(xd,yd) +% drawnow +% % axis equal +% plot(xd(model.nElements/4),yd(model.nElements/04), '*', 'Markersize' , 6) diff --git a/examples/vonKarman/model/S11.m b/examples/vonKarman/model/S11.m new file mode 100644 index 0000000..4715b9e --- /dev/null +++ b/examples/vonKarman/model/S11.m @@ -0,0 +1,80 @@ +function X = S11(model) +% This function assembles the quadratic nonlinearities for the curved beam +% example. + +E = model.E; +A = model.A; +dx = model.dx; +P = cell(1,6); + +P{1} =... +[ 0, 0, 0, 0, 0, 0;... + 0, -(3*A*E)/(5*dx^2), -(A*E)/(20*dx), 0, (3*A*E)/(5*dx^2), -(A*E)/(20*dx);... + 0, -(A*E)/(20*dx), -(A*E)/15, 0, (A*E)/(20*dx), (A*E)/60;... + 0, 0, 0, 0, 0, 0;... + 0, (3*A*E)/(5*dx^2), (A*E)/(20*dx), 0, -(3*A*E)/(5*dx^2), (A*E)/(20*dx);... + 0, -(A*E)/(20*dx), (A*E)/60, 0, (A*E)/(20*dx), -(A*E)/15]; + +P{2} =... +[ 0, -(3*A*E)/(5*dx^2), -(A*E)/(20*dx), 0, (3*A*E)/(5*dx^2), -(A*E)/(20*dx);... + -(3*A*E)/(5*dx^2), 0, 0, (3*A*E)/(5*dx^2), 0, 0;... + -(A*E)/(20*dx), 0, 0, (A*E)/(20*dx), 0, 0;... + 0, (3*A*E)/(5*dx^2), (A*E)/(20*dx), 0, -(3*A*E)/(5*dx^2), (A*E)/(20*dx);... + (3*A*E)/(5*dx^2), 0, 0, -(3*A*E)/(5*dx^2), 0, 0;... + -(A*E)/(20*dx), 0, 0, (A*E)/(20*dx), 0, 0]; + + P{3} = ... +[ 0, -(A*E)/(20*dx), -(A*E)/15, 0, (A*E)/(20*dx), (A*E)/60;... + -(A*E)/(20*dx), 0, 0, (A*E)/(20*dx), 0, 0;... + -(A*E)/15, 0, 0, (A*E)/15, 0, 0;... + 0, (A*E)/(20*dx), (A*E)/15, 0, -(A*E)/(20*dx), -(A*E)/60;... + (A*E)/(20*dx), 0, 0, -(A*E)/(20*dx), 0, 0;... + (A*E)/60, 0, 0, -(A*E)/60, 0, 0]; + +P{4} = ... +[ 0, 0, 0, 0, 0, 0;... + 0, (3*A*E)/(5*dx^2), (A*E)/(20*dx), 0, -(3*A*E)/(5*dx^2), (A*E)/(20*dx);... + 0, (A*E)/(20*dx), (A*E)/15, 0, -(A*E)/(20*dx), -(A*E)/60;... + 0, 0, 0, 0, 0, 0;... + 0, -(3*A*E)/(5*dx^2), -(A*E)/(20*dx), 0, (3*A*E)/(5*dx^2), -(A*E)/(20*dx);... + 0, (A*E)/(20*dx), -(A*E)/60, 0, -(A*E)/(20*dx), (A*E)/15]; + +P{5} = ... +[ 0, (3*A*E)/(5*dx^2), (A*E)/(20*dx), 0, -(3*A*E)/(5*dx^2), (A*E)/(20*dx);... + (3*A*E)/(5*dx^2), 0, 0, -(3*A*E)/(5*dx^2), 0, 0;... + (A*E)/(20*dx), 0, 0, -(A*E)/(20*dx), 0, 0;... + 0, -(3*A*E)/(5*dx^2), -(A*E)/(20*dx), 0, (3*A*E)/(5*dx^2), -(A*E)/(20*dx);... + -(3*A*E)/(5*dx^2), 0, 0, (3*A*E)/(5*dx^2), 0, 0;... + (A*E)/(20*dx), 0, 0, -(A*E)/(20*dx), 0, 0]; + +P{6} = ... +[ 0, -(A*E)/(20*dx), (A*E)/60, 0, (A*E)/(20*dx), -(A*E)/15;... + -(A*E)/(20*dx), 0, 0, (A*E)/(20*dx), 0, 0;... + (A*E)/60, 0, 0, -(A*E)/60, 0, 0;... + 0, (A*E)/(20*dx), -(A*E)/60, 0, -(A*E)/(20*dx), (A*E)/15;... + (A*E)/(20*dx), 0, 0, -(A*E)/(20*dx), 0, 0;... + -(A*E)/15, 0, 0, (A*E)/15, 0, 0]; + +S = cell(1,model.nDOF); +for i = 1:model.nDOF + S{i} = zeros(model.nDOF); +end +for i = 1:model.nElements + index = 3*(i-1)+1:3*(i+1); + Te = compute_rotation_matrix(model.nodes(i:i+1,:)); + for j = 1:6 + PPP = zeros(model.nDOF); + PPP(index,index) = Te.' * P{j} * Te; + ii = index(j); + S{ii} = S{ii}+ PPP; + end +end + +X = cell(1,length(model.freeDOFs)); +for j = 1:length(model.freeDOFs) + i = model.freeDOFs(j); + AJ = S{i}; + X{j} = AJ(model.freeDOFs,model.freeDOFs); +end + +end \ No newline at end of file diff --git a/examples/vonKarman/model/compute_rotation_matrix.m b/examples/vonKarman/model/compute_rotation_matrix.m new file mode 100644 index 0000000..b5b70be --- /dev/null +++ b/examples/vonKarman/model/compute_rotation_matrix.m @@ -0,0 +1,13 @@ +function T = compute_rotation_matrix(nodes) +T = speye(6,6); + +node1=nodes(1,:); +node2=nodes(2,:); +l=norm(node2-node1); +e_x = (node2-node1)/l; +R = [e_x(1) e_x(2); + -e_x(2) e_x(1)]; + +T(1:2,1:2) = R; +T(4:5,4:5) = R; +end \ No newline at end of file diff --git a/examples/vonKarman/model/sortVMsVK.m b/examples/vonKarman/model/sortVMsVK.m new file mode 100644 index 0000000..e36b37c --- /dev/null +++ b/examples/vonKarman/model/sortVMsVK.m @@ -0,0 +1,11 @@ +function [w,VN] = sortVMsVK(V,M,D) + +w = (diag(D)); +[w,pos] = sort(w); +VN = V(:,pos); + + +mu = diag(VN'*M*VN); +VN = VN*diag(1./sqrt(mu)); + +end diff --git a/modeselection/S11reduct.m b/modeselection/S11reduct.m new file mode 100644 index 0000000..3ccd40f --- /dev/null +++ b/modeselection/S11reduct.m @@ -0,0 +1,16 @@ +function S = S11reduct(SS,X) +% This function converts the quadratic nonlinearity coefficients from +% their general form to the desired form in modal coordinates + U = SS.U; + U2 = SS.U2; + B = size(U2); + S = cell(1,B(2)); + n = B(2)+SS.n; + + for j = 1:B(2) + S{j} = zeros(SS.n,SS.n); + for i = 1:n + S{j} = S{j} + U2(i,j) * U.' * X{i} * U; + end + end +end \ No newline at end of file diff --git a/modeselection/copyfromfull.m b/modeselection/copyfromfull.m new file mode 100644 index 0000000..22b2f69 --- /dev/null +++ b/modeselection/copyfromfull.m @@ -0,0 +1,15 @@ +function class = copyfromfull(class,full) +% copies properties of the second entry to the first entry +if ~isempty(full.f) + class.f = full.f; +end +if ~isempty(full.T) + class.T = full.T; + class.n_steps = full.n_steps; +end +if ~isempty(full.order) + class.order = full.order; +end + class.S = full.S; + class.DS = full.DS; +end \ No newline at end of file diff --git a/modeselection/modeselect.m b/modeselection/modeselect.m new file mode 100644 index 0000000..47ce78b --- /dev/null +++ b/modeselection/modeselect.m @@ -0,0 +1,101 @@ +function modes = modeselect(full,param,X) +%% Automated nonlinear mode selection procedure: described in Algorithm 1 of +% the following article +% G. Buza, S. Jain, G. Haller, Using Spectral Submanifolds for Optimal +% Mode Selection in Model Reduction, (2020) Preprint available on arXiv.org + + +try + modes = param.initialmodes; +catch + % initial set based on spectral quotient analysis + realpart = min(abs(real(-full.zeta.*full.omega0 + sqrt(full.zeta.^2-1).*full.omega0)),... + abs(real(-full.zeta.*full.omega0 - sqrt(full.zeta.^2-1).*full.omega0))); + [rpsorted,ind] = sort(realpart); + ratio = zeros(1,length(rpsorted)-1); + for i = 1:length(rpsorted)-1 + ratio(i) = rpsorted(i+1)/rpsorted(i); + end + cut = param.nmodes; + while cut > ceil(param.nmodes/3) % limit the number of modes selected via s.q.a. + [~,cut] = max(ratio); + if cut ~= 1 + ratio = ratio(1:cut-1); + end + end + modes = ind(1:cut).'; + + % refinement based on linear response to forcing (if given) + if ~isempty(full.f) + [x_linf,~] = full.LinearResponse(); + z = full.U.' * full.M * x_linf; + b = size(z); + NORM = zeros(1,b(1)); + for j = 1:b(1) + NORM(j) = norm(z(j,:)); + end + while length(modes) < ceil(param.nmodes*2/3) % limit the number of modes selected via linear response + if sum(NORM(modes)) > 0.9*sum(NORM) + break + end + modesC = setdiff(1:full.n,modes); + [~,ind2] = max(NORM(modesC)); + modes = [modes ind2]; + end + percent = 100*sum(NORM(modes))/sum(NORM); + disp(['The initial mode selection recovers ' num2str(percent) '% of the linear response.']) + else + disp('Since no forcing was provided for the SSR class, it is not guaranteed that the mode selection recovers the linear response.') + end +end +disp('Initial linear selection:') +disp(modes) + +% nonlinear selection +itnumber = 0; +while length(modes) < param.nmodes + red = SSR(full.M,full.C,full.K,modes); + red = copyfromfull(red,full); % copies general props + switch param.example % obtaining the nonlinear coefficients is problem dependent + case 'simple' % heavy on memory + R = S11reduct(red,X); + case 'plate' % memory preserving approach + V = zeros(param.fulldof,red.n); + VC = zeros(param.fulldof,length(red.M)-red.n); + V(param.freedofs,:) = red.U; + VC(param.freedofs,:) = red.U2; + [K2,~] = AssembleTensors(param.model,V,VC); % computes quadratic coeffs R_i + R = cell(1,size(red.U2,2)); + for i = 1:size(red.U2,2) + R{i} = K2(i,:,:); + end + end + W = red.SSM2(R); % computing SSM + A = [zeros(red.n),eye(red.n);-red.K1,-red.C1]; + [~ ,scaldir] = scalg(A,W); % computing directional curvatures + [scnorm,scindex] = sort(abs(scaldir)); + Jmodes = setdiff(1:full.n,modes); + scindex = fliplr(scindex); + crit = sum(scnorm) * (1-param.curvtolerance); + Z = 0; + numrec = 0; + while Z < crit + Z = Z + max(scnorm); + scnorm = setdiff(scnorm,max(scnorm)); + numrec = numrec + 1; + end + recomm = Jmodes(scindex(1:numrec)); % recommnended modes + check = param.nmodes - (length(modes) + length(recomm)); % check if it overshoots the number of allowed modes + if check < 0 + modes = sort([modes recomm(1:(param.nmodes - length(modes)))]); + else + modes = sort([modes recomm]); + end + disp(modes) + itnumber = itnumber + 1; + if ~param.type + break + end +end +disp(['The nonlinear mode selection procedure was performed ' num2str(itnumber) ' time(s).']) +end \ No newline at end of file diff --git a/modeselection/scalg.m b/modeselection/scalg.m new file mode 100644 index 0000000..919105c --- /dev/null +++ b/modeselection/scalg.m @@ -0,0 +1,24 @@ +function [scal,scaldir] = scalg(A,W) +% This function computes directional scalar curvatures. + +m = length(W{1}); + +SUM = 0; +scaldir = zeros(1,length(W)); +for l = 1:length(W) + Wl = W{l}; + sumk = 0; + for i = 1:m + for j = 1:m + a = 4 * ( Wl(i,i) * Wl(j,j) - (Wl(i,j))^2 +... + 4* (Wl(i,:)*A(:,i))*(Wl(j,:)*A(:,j))-... + (Wl(i,:)*A(:,j) + Wl(j,:)*A(:,i))^2); + SUM = SUM + a; + sumk = sumk+a; + end + end + scaldir(l) = sumk; +end + +scal = SUM/(m^2); +end \ No newline at end of file