Skip to content

[WIP] Energy preserving #67

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 19 commits into
base: master
Choose a base branch
from
Open

[WIP] Energy preserving #67

Changes from 17 commits
Commits
Show all changes
19 commits
Select commit Hold shift + click to select a range
6ab5752
Energy-preserving non linear conditions
Sondar74 Sep 19, 2023
fbd9bdd
Nonlinear constraints
Sondar74 Sep 24, 2023
c1fe866
Conditions for PEP 7
Sondar74 Nov 21, 2023
53a9a97
Eliminate duplication
Sondar74 Nov 22, 2023
7a8eaa1
Fix PEP 7th order
Sondar74 Nov 23, 2023
48eebed
Nonlinear constraints PEP order 7
Sondar74 Dec 5, 2023
7b06e40
PEP7
Sondar74 Dec 25, 2023
e252998
instructions
Sondar74 Dec 25, 2023
b2f1fa3
restore
Sondar74 Dec 25, 2023
e28e460
restore
Sondar74 Dec 25, 2023
ab74270
Generate q=7 constraints
Sondar74 Jan 16, 2024
bb658bf
Project and Manifest
Sondar74 Mar 14, 2024
0f84f6f
Update RK-coeff-opt/PEP_Generator/README.txt
Sondar74 Mar 14, 2024
8e48274
closing PR
Sondar74 Apr 15, 2024
304fd03
Merge branch 'EnergyPreserving' of github.com:Sondar74/RK-Opt into En…
Sondar74 Apr 15, 2024
7c2bba3
Update RK-coeff-opt/nonlinear_constraints.m
Sondar74 Apr 16, 2024
e57a2c5
Error message
Sondar74 Apr 18, 2024
c1fd47c
Move message to the end
Sondar74 Apr 19, 2024
ec04e9c
Update conditions 6 and 7
Sondar74 May 7, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
63 changes: 63 additions & 0 deletions RK-coeff-opt/nonlinear_constraints.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -120,3 +120,66 @@
res = polyval(poly_coef, constrain_emb_stability);
con = [con, res .* conj(res) - 1];
end

q=1;
% Throw error message in case q>7
msg = 'The energy-preserving conditions for orders q>7 have not been implemented.';
% PEP conditions
if q >= 8
error(msg);
end
if q>=1
coneq(end+1) = sum(b)-1;
end
if q>=2
coneq(end+1) = c'*b-1/2;
end
if q >= 3
coneq(end+1) = c'.^2*b-1/3;
end
if q >= 4
coneq(end+1) = b'*A^2*c-b'*A*c+1/8;
coneq(end+1) = c'.^3*b-1/4;
coneq(end+1) = ((b'.*c')*A*c)-(b'*A*c.^2)/2-1/12;
end
if q >= 5
coneq(end+1) = (b'*A^2*c.^2)/2 + (b.*c)'*A^2*c-b'*A^2*c-(b'*A*c.^2)/2 + (b'*A*c)/2 - 1/24;
coneq(end+1) = c'.^4*b-1/5;
coneq(end+1) = 2*(b'*A*diag(c)*A*c)-(b'*(A*c).^2) - (b'*A^2*c)-(b'*A*c.^2)+(b'*A*c)-1/24;
coneq(end+1) = (b.*c.^2)'*A*c-(b'*A*c.^3)/3 - 1/12;
end
if q >= 6
coneq(end+1) = -(b'*A*c)^2 + (b'*A^2*c)- 2 *(b'*A^3*c)+ 2* (b'*A^4*c);
coneq(end+1) = 6* (3* (b'*A*c)+ ((b'.*c')*A*c)+ 3* ( (b'*A^2*c.^2)+ 4* (b'*diag(c)*A^3*c))) - 2 - 3* (b'*A*c.^2)- 24* (b'*A^2*c)- 36* ((b.*c)'*A^2*c)- 36* (b'*A^3*c.^2);
coneq(end+1) = 1 - 12* (b'*A*c)+ 24* ((b'.*c')*A*c)- 36* ((b.*c.^2)'*A*c)+ 24* (b'*A*c.^2)- 12* (b'*A*c.^3)+ 24* (b'*A^2*c)- 72* ((b.*c)'*A^2*c)+ 72* (b'*diag(c).^2*A^2*c)- 36* (b'*A^2*c.^2)+ 24* (b'*A^2*c.^3);
coneq(end+1) = 12 *(2* ((b'.*c')*A*c)+ 6* ((b.*c.^2)'*A*c)+ 6* (b'*(A*c).^2)+ 5* (b'*A*c.^2)+ 2* (b'*A^2*c)+ 12* b'*A*diag(c).^2*A*c) - 5 - 24* (b'*A*c)- 144* (b'*diag(c)*(A*c).^2)- 72* (b'*A*c.^3)- 144* (b'*A*diag(c)*A*c);
coneq(end+1) = 1 + 60* ((b'.*c')*A*c)+ 120* (b'*diag(c).^3*A*c)+ 60* (b'*A*c.^3)- 30* (6* ((b.*c.^2)'*A*c)+ (b'*A*c.^2)+ (b'*A*c.^4));
coneq(end+1) = 1 + 12* ((b'.*c')*A*c)+ 36* (b'*diag(c).^2*A*c.^2)+ 12* (b'*A*c.^3)-6* (6* ((b.*c.^2)'*A*c)+ (b'*A*c.^2)+ 4* (b'*diag(c)*A*c.^3));
coneq(end+1) = -1 + 6 * (b'*A*c)+ 6* ((b'.*c')*A*c)- 3* (b'*A*c.^2)+ 12* (b'*A^2*c)- 36* ((b.*c)'*A^2*c)- 18* (b'*A^2*c.^2)+ 36* (b'*diag(c)*A^2*c.^2);
coneq(end+1) = -(5/144) - (b'*A*c)/6 + (2* ((b'.*c')*A*c))/3 - (b'*(A*c).^2)/2 + (b'*A*c.^2)/6 + (b'*A^2*c)/6 - (b'*A*diag(c)*A*c)+ (b'*A*(A*c).^2);
coneq(end+1) = c'.^5*b-1/6;
coneq(end+1) = 8 *( (b'*A*c)+ 3* ( (b'*(A*c).^2)+ (b'*A^2*c.^2)+ 2* (b'*A*diag(c)*A^2*c))) - 1 - 12* (b'*A*c.^2)- 8 * (b'*A^2*c)- 48* (b'*((A^2*c).*(A*c)))-48* (b'*A^2*diag(c)*A*c);
coneq(end+1) = 2* (6* ((b'.*c')*A*c)+ 12* (b'*(A*c).^2)+ 9* (b'*A*c.^2)+ 8* (b'*A^2*c)+ 24* (b'*diag(c)*A*diag(c)*A*c)+ 12* (b'*A*diag(c)*A*c.^2)) -1 - 4* (b'*A*c)- 24* ((b.*c)'*A*c.^2)- 24* (b'*((A*c.^2).*(A*c)))- 24* ((b.*c)'*A^2*c)- 48* (b'*A*diag(c)*A*c)- 12* (b'*A^2*c.^2);
end
if q >= 7
coneq(end+1) = -(b'*A^4*c) + (b'*((A*(A*(A*(A*c)))).*c)) + (b'*((A*(A*(A*(A*c.^2))))))/2 + (b'*A^3*c)/2 - (b'*diag(c)*A^3*c)/2 - (b'*A^3*c.^2)/4 - (b'*A^2*c)/12 + ((b.*c)'*A^2*c)/12 + (b'*A^2*c.^2)/24 - ((b'*A*c)/6 - 1/36)*((b'*A*c) - 1/6) - ((b'*A*c)/3 - 1/18)*((b'*A*c) - 1/6) + ((b'*A*c)/2 - 1/12)*((b'*A*c) - 1/6) - ((b'*A*c)/2 - 1/12)*(-(b'*A*c)/2 + ((b'.*c')*A*c) - 1/24) - ((b'*A*c) - 1/6)*(-(b'*A*c)/2 + c'.^3*b/2 + 1/24)/2;
coneq(end+1) = -(b'*A^4*c)/2 + (b'*((A*(A*diag(c)*(A*(A*c)))))) + (b'*A^3*c)/2 - (b'*A*diag(c)*A^2*c)/2 - (b'*A^2*diag(c)*A*c)/2 + (b'*((A*(A*c)).*(A*(A*c))))/2 - (b'*((A^2*c).*(A*c)))/2 - (b'*A^2*c)/4 + ((b.*c)'*A^2*c)/6 + (b'*A*diag(c)*A*c)/3 + (b'*A^2*c.^2)/12 + (b'*(A*c).^2)/12 + (b'*A*c)/12 - ((b'.*c')*A*c)/12 - (b'*A*c.^2)/12 - ((b'*A*c)/6 - 1/36)*((b'*A*c) - 1/6) - 5*((b'*A*c)/3 - 1/18)*((b'*A*c) - 1/6)/2 + 3*((b'*A*c)/2 - 1/12)*((b'*A*c) - 1/6)/2 - 2*((b'*A*c)/2 - 1/12)*((b'*A^2*c) - (b'*A*c) + 1/8);
coneq(end+1) = (b'*diag(c)*A^3*c)/2 - (b'*((A*(A*(A*c))).*c.^2))/2 - (b'*A^3*c.^2)/4 + (b'*((A*(A*(A*c.^3)))))/6 - (b'*A^2*c)/12 - ((b.*c)'*A^2*c)/6 + (b'*diag(c).^2*A^2*c)/4 + (b'*A^2*c.^2)/6 - (b'*A^2*c.^3)/12 + (b'*A*c)/24 - ((b.*c.^2)'*A*c)/24 - (b'*A*c.^2)/24 + (b'*A*c.^3)/72 + ((b'*A*c)/2 - 1/12)*((c'.^3*b) - 1/4)/2 - ((b'*A*c) - 1/6)*((c'.^3*b)/2 - 1/8)/6;
coneq(end+1) = -((b.*c.^2)'*A*c)/12 + (b'*diag(c).^3*A*c)/12 + (b'*diag(c).^2*A*c.^2)/8 - (b'*((A*c.^2).*c.^3))/12 + (b'*A*c.^3)/36 - (b'*diag(c)*A*c.^3)/12 - (b'*A*c.^4)/48 + (b'*((A*c.^4).*c))/24 + c'.^3*b/72 - 5*(c'.^4*b)/288;
coneq(end+1) = b'*((A*(A*((A*c).*(A*c)))))/2 - (b'*((A*((A*c).*(A*(A*c)))))) + (b'*A*diag(c)*A^2*c)/2 - (b'*A^2*diag(c)*A*c)/2 + (b'*A*(A*c).^2)/4 + (b'*((A^2*c).*(A*c)))/2 - ((b.*c)'*A^2*c)/3 + (b'*A*diag(c)*A*c)/12 + (b'*A^2*c.^2)/12 - 7*(b'*(A*c).^2)/24 + (b'*A*c)/24 + ((b'.*c')*A*c)/6 - (b'*A*c.^2)/12 + ((b'*A*c)/6 - 1/36)*((b'*A*c) - 1/6) + ((b'*A*c)/3 - 1/18)*((b'*A*c) - 1/6) - ((b'*A*c)/2 - 1/12)*((b'*A*c) - 1/6) + ((b'*A*c)/2 - 1/12)*((b'*A^2*c) - (b'*A*c) + 1/8) - 1/144;
coneq(end+1) = (b'*A^3*c)/6 - (b'*A*diag(c)*A^2*c)/2 - (b'*A^2*diag(c)*A*c)/2 + (b'*((A*diag(c)*(A*diag(c)*(A*c))))) - (b'*((A^2*c).*(A*c)))/2 - (b'*A^2*c)/6 + ((b.*c)'*A^2*c)/2 + (b'*((A*diag(c)*(A*c)).*(A*c))) + 2*(b'*A*diag(c)*A*c)/3 - (b'*diag(c)*A*diag(c)*A*c) + (b'*A^2*c.^2)/6 - (b'*A*diag(c)*A*c.^2)/2 + (b'*(A*c).^2)/4 - (b'*A*c)/180 - ((b'.*c')*A*c)/4 - (b'*((A*c.^2).*(A*c)))/2 - (b'*A*c.^2)/6 + ((b.*c)'*A*c.^2)/2 - 2*((b'*A*c)/6 - 1/36)*((b'*A*c) - 1/6) + ((b'*A*c)/2 - 1/12)*((b'*A*c) - 1/6) - 2*((b'*A*c) - 1/6)*(-(b'*A*c)/4 + ((b'.*c')*A*c)/2 - 1/48) + 149/15120;
coneq(end+1) = ((b.*c)'*A^2*c)/12 - (b'*diag(c).^2*A^2*c)/4 + (b'*((A*(A*c)).*c.^3))/6 + (b'*A^2*c.^2)/24 - (b'*A^2*c.^3)/12 + (b'*((A*(A*c.^4))))/24 + ((b.*c.^2)'*A*c)/12 - (b'*diag(c).^3*A*c)/12 - (b'*A*c.^2)/24 + (b'*A*c.^3)/18 - (b'*A*c.^4)/48 - (c'.^3*b)/72 + 5*(c'.^4*b)/288 - ((b'*A*c)/2 - 1/12)*((c'.^3*b) - 1/4)/6;
coneq(end+1) = ((b.*c)'*A^2*c)/4 - (b'*diag(c).^2*A^2*c)/4 + (b'*A^2*c.^2)/8 - (b'*diag(c)*A^2*c.^2)/2 + (b'*((A*(A*c.^2)).*c.^2))/4 - (b'*A^2*c.^3)/12 + (b'*((A*(A*c.^3)).*c))/6 - (b'*A*c)/24 + ((b.*c.^2)'*A*c)/24 + (b'*A*c.^2)/24 - (b'*A*c.^3)/72;
coneq(end+1) = -3*(b'*A*(A*c).^2)/4 + (b'*((A*((A*c).*(A*c))).*c))/2 + ((b.*c)'*A^2*c)/12 + 5*(b'*A*diag(c)*A*c)/12 - (b'*diag(c)*A*diag(c)*A*c)/2 + (b'*((A*(c.^2.*(A*(A*c))))))/2 + (b'*A^2*c.^2)/24 - (b'*A*diag(c)*A*c.^2)/4 + 7*(b'*(A*c).^2)/24 - (b'*A*c)/24 - ((b'.*c')*A*c)/6 - (b'*((A*c.^2).*(A*c)))/4 - (b'*A*c.^2)/24 + ((b.*c)'*A*c.^2)/4 + 1/144;
coneq(end+1) = ((b.*c)'*A^2*c)/6 + (b'*A*diag(c)*A*c)/6 - (b'*diag(c)*A*diag(c)*A*c)/2 + (b'*A^2*c.^2)/12 - (b'*diag(c)*A^2*c.^2)/4 - (b'*A*diag(c)*A*c.^2)/4 + (b'*((A*(c.*(A*c.^2))).*c))/2 - (b'*(A*c).^2)/12 - (b'*A*c)/12 + ((b'.*c')*A*c)/12 + (b'*((A*c.^2).*(A*c)))/4 - (b'*((A*c.^2).*(A*c.^2)))/8 + (b'*A*c.^2)/24;
coneq(end+1) = -(b'*A*(A*c).^2)/4 + (b'*((A*(c.*(A*c).*(A*c)))))/2 + (b'*A*diag(c)*A*c)/4 - (b'*A*diag(c).^2*A*c)/2 - (b'*((A*c).*(A*c).*(A*c)))/6 + (b'*(A*c).^2)/8 - ((b'.*c')*A*c)/4 + ((b.*c.^2)'*A*c)/4 + (b'*A*c.^3)/12 - (c'.^3*b)/12 + 1/48;
coneq(end+1) = (b'*A*diag(c)*A*c)/12 - (b'*A*diag(c).^2*A*c)/4 + (b'*((A*(c.^3.*(A*c)))))/6 - (b'*(A*c).^2)/24 + (b'*diag(c)*(A*c).^2)/4 - (b'*((A*c).*(A*c).*c.^2))/4 - ((b'.*c')*A*c)/12 - ((b.*c.^2)'*A*c)/24 + (b'*diag(c).^3*A*c)/6 + 7*(b'*A*c.^3)/72 - (b'*A*c.^4)/12 - 7*(c'.^3*b)/72 + (c'.^4*b)/72 + 1/48;
coneq(end+1) = -((b.*c.^2)'*A*c)/24 + (b'*diag(c).^3*A*c)/12 - (b'*((A*c).*c.^4))/24 + (b'*A*c.^3)/72 - (b'*A*c.^4)/48 + (b'*((A*c.^5)))/120 - (c'.^3*b)/72 - 5*(c'.^4*b)/288 + (c'.^5*b)/60 + 1/240;
coneq(end+1) = -(b'*A^4*c) + (b'*((A*(c.*(A*(A*(A*c))))))) + (b'*((A*(A*(A*(c.*(A*c))))))) - (b'*((A*(A*(A*c))).*(A*c))) + (b'*A^3*c) - (b'*A*diag(c)*A^2*c)/2 - (b'*A^2*diag(c)*A*c)/2 - (b'*A^3*c.^2)/2 + (b'*((A^2*c).*(A*c)))/2 - (b'*A^2*c)/2 + ((b.*c)'*A^2*c)/6 + (b'*A*diag(c)*A*c)/6 + (b'*A^2*c.^2)/3 - (b'*(A*c).^2)/12 + (b'*A*c)/8 - ((b'.*c')*A*c)/12 - (b'*A*c.^2)/12 - ((b'*A*c)/6 - 1/36)*((b'*A*c) - 1/6) - ((b'*A*c)/3 - 1/18)*((b'*A*c) - 1/6) + ((b'*A*c)/2 - 1/12)*((b'*A*c) - 1/6) - ((b'*A*c)/2 - 1/12)*(-(b'*A*c) + (b'*A*c.^2) + 1/12) + ((b'*A*c)/2 - 1/12)*(-(b'*A*c)/2 + ((b'.*c')*A*c) - 1/24) - ((b'*A*c) - 1/6)*(-(b'*A*c)/4 + ((b'.*c')*A*c)/2 - 1/48);
coneq(end+1) = (b'*diag(c)*A^3*c)/2 + (b'*A*diag(c)*A^2*c)/2 - (b'*((A*(c.*(A*(A*c)))).*c)) - (b'*A^2*diag(c)*A*c)/2 - (b'*A^3*c.^2)/4 + (b'*((A*(A*(c.*(A*c.^2))))))/2 - (b'*((A^2*c).*(A*c)))/2 + (b'*((A*(A*c)).*(A*c.^2)))/2 - (b'*A^2*c)/4 + (b'*A*diag(c)*A*c)/6 + (b'*diag(c)*A*diag(c)*A*c)/2 + (b'*A^2*c.^2)/4 - (b'*A*diag(c)*A*c.^2)/4 + (b'*(A*c).^2)/6 + (b'*A*c)/6 - ((b'.*c')*A*c)/6 - (b'*((A*c.^2).*(A*c)))/4 - (b'*A*c.^2)/8 - ((b'*A*c)/2 - 1/12)*(-(b'*A*c) + (b'*A*c.^2) + 1/12)/2 + ((b'*A*c)/2 - 1/12)*(-(b'*A*c)/2 + ((b'.*c')*A*c) - 1/24);
coneq(end+1) = (b'*diag(c)*A^3*c)/2 - (b'*A*diag(c)*A^2*c)/2 + (b'*A^2*diag(c)*A*c)/2 - (b'*((A*(A*(c.*(A*c)))).*c)) - (b'*A^3*c.^2)/4 + (b'*((A*(c.*(A*(A*c.^2))))))/2 + (b'*((A^2*c).*(A*c)))/2 - (b'*A^2*c)/4 - ((b.*c)'*A^2*c)/12 + (b'*A*diag(c)*A*c)/6 - (b'*((A*(A*c.^2)).*(A*c)))/2 - (b'*A^2*c.^2)/24 + (b'*diag(c)*A^2*c.^2)/2 - (b'*(A*c).^2)/12 + (b'*A*c)/8 - ((b'.*c')*A*c)/12 - (b'*A*c.^2)/12;
coneq(end+1) = (b'*A^3*c)/6 - (b'*A*diag(c)*A^2*c)/2 + (b'*((A*(c.^2.*(A*(A*c))))))/2 - (b'*A^2*diag(c)*A*c)/2 + (b'*((A*(A*(c.^2.*(A*c))))))/2 + (b'*((A^2*c).*(A*c)))/2 - (b'*((A*(A*c)).*(A*c).*c)) - (b'*A^2*c)/6 + ((b.*c)'*A^2*c)/12 + (b'*diag(c).^2*A^2*c)/4 + (b'*A*diag(c)*A*c)/2 - (b'*A*diag(c).^2*A*c)/2 + 5*(b'*A^2*c.^2)/24 - (b'*A^2*c.^3)/4 - (b'*(A*c).^2)/6 + (b'*diag(c)*(A*c).^2)/2 + 3*(b'*A*c)/40 - ((b'.*c')*A*c)/4 - ((b.*c.^2)'*A*c)/24 - 5*(b'*A*c.^2)/24 + 5*(b'*A*c.^3)/24 - (c'.^3*b)/12 - ((b'*A*c)/2 - 1/12)*(-(b'*A*c) + (b'*A*c.^2) + 1/12)/2 + ((b'*A*c)/2 - 1/12)*(-(b'*A*c)/2 + ((b'.*c')*A*c) - 1/24) + 59/1680;
coneq(end+1) = ((b.*c)'*A^2*c)/6 - (b'*diag(c).^2*A^2*c)/4 + (b'*A*diag(c)*A*c)/6 - (b'*diag(c)*A*diag(c)*A*c)/2 + (b'*((A*(c.*(A*c))).*c.^2))/2 + (b'*A^2*c.^2)/12 - (b'*A*diag(c)*A*c.^2)/4 - (b'*A^2*c.^3)/12 + (b'*((A*(c.*(A*c.^3)))))/6 - (b'*(A*c).^2)/12 - ((b'.*c')*A*c)/12 + ((b.*c.^2)'*A*c)/8 + (b'*((A*c.^2).*(A*c)))/4 - (b'*A*c.^2)/12 + ((b.*c)'*A*c.^2)/4 - (b'*diag(c).^2*A*c.^2)/4 - (b'*((A*c.^3).*(A*c)))/6 + (b'*A*c.^3)/24;
coneq(end+1) = ((b.*c)'*A^2*c)/12 + (b'*A*diag(c)*A*c)/3 - (b'*diag(c)*A*diag(c)*A*c)/2 - (b'*A*diag(c).^2*A*c)/2 + (b'*((A*(c.^2.*(A*c))).*c))/2 + (b'*A^2*c.^2)/24 - (b'*A*diag(c)*A*c.^2)/4 + (b'*((A*(c.^2.*(A*c.^2)))))/4 - (b'*(A*c).^2)/6 + (b'*diag(c)*(A*c).^2)/2 - (b'*A*c)/24 - ((b'.*c')*A*c)/12 - ((b.*c.^2)'*A*c)/8 + (b'*((A*c.^2).*(A*c)))/4 - (b'*((A*c.^2).*(A*c).*c))/2 - (b'*A*c.^2)/24 + ((b.*c)'*A*c.^2)/4 + (b'*diag(c).^2*A*c.^2)/8 + (b'*A*c.^3)/8 - (b'*diag(c)*A*c.^3)/4 + 1/144;
coneq(end+1) = (c'.^4*b)/288 - (c'.^5*b)/240 + b'*(c.^6)/720 - 1/5040;
end