-
Notifications
You must be signed in to change notification settings - Fork 4
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
* Modeselection * minor changes #1 * ensuring sparse/full eigenvalue computation in the reduced/full case. * Commentary Co-authored-by: gergelybuza <[email protected]>
- Loading branch information
1 parent
12dc06e
commit f6704e3
Showing
30 changed files
with
1,104 additions
and
15 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -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 |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -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 |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -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 |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -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 |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Oops, something went wrong.