clc; clear;
echo on
pvar s theta;
pde_struct PDE;
a = 0; b = 1;
PDE.x{1}.vars = s;
PDE.x{1}.dom = [a,b];
PDE.x{1}.diff = 2;
PDE.x{1}.term{1}.D = 1;
PDE.BC{1}.term{1}.loc = a;
PDE.BC{2}.term{1}.loc = b;
PDE = initialize(PDE);
PIE = convert(PDE,'pie');
H2 = PIE.T;
H1 = PIE.A;
dpvar gam;
vars = [H2.var1; H2.var2];
prob = sosprogram(vars,gam);
prob = sossetobj(prob, gam);
opts.psatz = 1;
prob = lpi_ineq(prob,-(H2'*H2-gam*H1'*H1),opts);
prob = sossolve(prob);
poincare_constant = sqrt(double(sosgetsol(prob,gam)));
echo off
fprintf(['\n If successful, ',num2str(poincare_constant),' is an upper bound on Poincare''s constant for this problem.\n'])
fprintf([' An optimal value of Poincare''s constant on domain [0,1] is known to be 1/pi=',num2str(1/(pi)),'.\n']);
%%%%%%%%%%%%%%%%%% Start Code Snippet %%%%%%%%%%%%%%%%%%
%%%%% Declare the PDE
% % Initialize the PDE structure and spatial variable s in [a,b]
pvar s theta;
pde_struct PDE;
a = 0; b = 1;
% % Declare the state variables x(t,s)
PDE.x{1}.vars = s;
PDE.x{1}.dom = [a,b];
PDE.x{1}.diff = 2; % Let x be second order differentiable wrt s.
% % Declare the PDE \dot{x}(t,s) = \partial_{s} x(t,s)
PDE.x{1}.term{1}.D = 1; % Order of the derivative wrt s
% % Declare the boundary conditions x(t,a) = x(t,b) = 0
PDE.BC{1}.term{1}.loc = a; % Evaluate x at s=a
PDE.BC{2}.term{1}.loc = b; % Evaluate x at s=b
% % Initialize the system
PDE = initialize(PDE);
Encountered 1 state component:
x(t,s), of size 1, differentiable up to order (2) in variables (s);
Encountered 2 boundary conditions:
F₁(t) = 0, of size 1;
F₂(t) = 0, of size 1;
%%%%% Convert the PDE to a PIE
PIE = convert(PDE,'pie');
--- Reordering the state components to allow for representation as PIE ---
--- Converting ODE-PDE to PIE ---
H2 = PIE.T; % H2 x_{ss} = x
H1 = PIE.A; % H1 x_{ss} = x_{s}
%%%%% Solve the LPI < H2 x_ss, H2 x_ss> - gam <H1 x_ss, H1 x_ss> <=0
%%%%% where (H2 x_ss) = x and (H1 x_ss) = x_s
% % First, define dpvar gam and set up an optimization problem
dpvar gam;
vars = [H2.var1; H2.var2];
prob = sosprogram(vars,gam);
% % Set gam as objective function to minimize
prob = sossetobj(prob, gam);
% % Set up the constraint H2'*H2-gam H1'*H1<=0
opts.psatz = 1; % Add psatz term to allow H2'*H2 > gam H1'*H1 outside of [a,b]
prob = lpi_ineq(prob,-(H2'*H2-gam*H1'*H1),opts);
% Solve and retrieve the solution
prob = sossolve(prob);
Size: 339 69
The coefficient matrix is not full row rank, numerical problems may occur.
SeDuMi 1.3 by AdvOL, 2005-2008 and Jos F. Sturm, 1998-2003.
Alg = 2: xz-corrector, Adaptive Step-Differentiation, theta = 0.250, beta = 0.500
Put 1 free variables in a quadratic cone
eqs m = 69, order n = 29, dim = 341, blocks = 4
nnz(A) = 1023 + 0, nnz(ADA) = 4709, nnz(L) = 2389
it : b*y gap delta rate t/tP* t/tD* feas cg cg prec
0 : 1.01E+00 0.000
1 : 2.94E-01 2.66E-01 0.000 0.2632 0.9000 0.9000 0.51 1 1 2.5E+00
2 : 4.11E-01 8.60E-02 0.000 0.3231 0.9000 0.9000 0.25 1 1 1.2E+00
3 : 3.63E-01 2.23E-02 0.000 0.2600 0.9000 0.9000 0.64 1 1 3.3E-01
4 : 2.25E-01 6.40E-03 0.000 0.2862 0.9000 0.9000 1.23 1 1 8.2E-02
5 : 1.61E-01 1.73E-03 0.000 0.2704 0.9000 0.9000 1.28 1 1 1.9E-02
6 : 1.49E-01 4.89E-04 0.000 0.2825 0.9000 0.9000 1.14 1 1 5.1E-03
7 : 1.57E-01 1.72E-04 0.000 0.3528 0.9000 0.9000 0.81 1 1 2.1E-03
8 : 1.64E-01 7.12E-05 0.000 0.4132 0.9000 0.9000 0.64 1 1 1.1E-03
9 : 1.68E-01 2.78E-05 0.000 0.3905 0.9000 0.9000 0.69 1 1 4.9E-04
10 : 1.71E-01 1.37E-05 0.000 0.4917 0.9000 0.9000 0.58 1 1 3.1E-04
11 : 1.73E-01 5.79E-06 0.000 0.4235 0.9000 0.9000 0.61 1 1 1.6E-04
12 : 1.75E-01 3.09E-06 0.000 0.5333 0.9000 0.9000 0.48 1 1 1.1E-04
13 : 1.76E-01 1.26E-06 0.000 0.4077 0.9000 0.9000 0.58 1 1 5.5E-05
14 : 1.78E-01 6.30E-07 0.000 0.5005 0.9000 0.9000 0.43 1 1 3.8E-05
15 : 1.79E-01 2.52E-07 0.000 0.3991 0.9000 0.9000 0.54 1 2 1.9E-05
16 : 1.79E-01 1.22E-07 0.000 0.4852 0.9000 0.9000 0.41 1 2 1.3E-05
17 : 1.80E-01 4.86E-08 0.000 0.3978 0.9000 0.9000 0.54 2 2 6.4E-06
18 : 1.80E-01 1.08E-08 0.000 0.2230 0.6845 0.9000 0.41 2 2 4.3E-06
19 : 1.81E-01 3.65E-09 0.000 0.3373 0.9000 0.9135 0.53 2 2 1.9E-06
20 : 1.81E-01 1.62E-09 0.000 0.4425 0.9081 0.9000 0.51 3 3 1.1E-06
21 : 1.81E-01 9.01E-10 0.000 0.5579 0.6076 0.9000 0.44 5 4 8.4E-07
22 : 1.81E-01 3.24E-10 0.000 0.3593 0.9000 0.9175 0.57 5 4 4.0E-07
23 : 1.81E-01 1.44E-10 0.000 0.4447 0.9143 0.9000 0.46 5 5 2.4E-07
24 : 1.82E-01 8.36E-11 0.000 0.5807 0.5466 0.9000 0.37 5 6 1.9E-07
25 : 1.82E-01 2.83E-11 0.000 0.3386 0.9000 0.9216 0.46 7 6 9.3E-08
26 : 1.82E-01 1.34E-11 0.000 0.4727 0.9100 0.9000 0.42 7 8 6.0E-08
27 : 1.82E-01 7.63E-12 0.000 0.5702 0.5890 0.9000 0.32 10 10 4.8E-08
28 : 1.82E-01 3.05E-12 0.000 0.3991 0.9000 0.9179 0.48 12 10 2.5E-08
29 : 1.82E-01 1.44E-12 0.000 0.4720 0.8628 0.9000 0.39 18 21 1.7E-08
30 : 1.82E-01 8.03E-13 0.000 0.5584 0.9000 0.9000 0.39 15 14 1.3E-08
31 : 1.82E-01 5.58E-13 0.000 0.6943 0.9000 0.9489 0.45 15 15 1.0E-08
32 : 1.82E-01 2.30E-13 0.000 0.4131 0.9000 0.9432 0.52 24 25 5.1E-09
iter seconds digits c*x b*y
32 0.6 3.6 1.8203777777e-01 1.8208352246e-01
|Ax-b| = 1.1e-08, [Ay-c]_+ = 3.4E-11, |x|= 5.2e+03, |y|= 6.8e+04
Detailed timing (sec)
Pre IPM Post
4.700E-02 2.280E-01 1.599E-02
Max-norms: ||b||=5.000000e-01, ||c|| = 1,
Cholesky |add|=3, |skip| = 35, ||L.L|| = 1.64194e+07.
Residual norm: 1.1069e-08
iter: 32
feasratio: 0.5198
pinf: 0
dinf: 0
numerr: 0
timing: [0.0470 0.2280 0.0160]
wallsec: 0.2910
cpusec: 0.6719
poincare_constant = sqrt(double(sosgetsol(prob,gam)));
%%%%%%%%%%%%%%%%%% End Code Snippet %%%%%%%%%%%%%%%%%%
echo off
If successful, 0.42666 is an upper bound on Poincare's constant for this problem.
An optimal value of Poincare's constant on domain [0,1] is known to be 1/pi=0.31831.