designingpolynomials.m

% Data 
x0 = -1
x1 = 0;
x2 = 1;
y0 = 0;
y1 = 1;
y2 = 0;
n = 9;

% Solve conditions
x = sdpvar(1);
[p,a,v] = polynomial(x,n);
Model = [replace(p,x,x0)==y0,
         replace(p,x,x1)>=y1,
         replace(p,x,x2)==y2];
optimize(Model)

xv = linspace(x0,x2,100);
yv = polyval(fliplr(value(a')),xv);
clf
plot(xv,yv)

% Add derivative conditions
dp = jacobian(p,x);
dp2 = hessian(p,x);
Model = [replace(p,x,x0)==y0,
         replace(p,x,x1)>=y1, 
         replace(p,x,x2)==y2,
         replace(dp,x,x0)==0,
         replace(dp,x,x2)==0,
         replace(dp2,x,x1)>=0];
optimize(Model)
yv = polyval(fliplr(value(a')),xv);
hold on
plot(xv,yv)

%% Add integral objective
Model = [replace(p,x,x0)==y0,
         replace(p,x,x1)>=y1, 
         replace(p,x,x2)==y2,
         replace(dp,x,x0)==0,
         replace(dp,x,x2)==0,
         replace(dp2,x,x1)>=0];
optimize(Model, int(p^2,x,-1,1));
yv = polyval(fliplr(value(a')),xv);
hold on
plot(xv,yv)

%% Gridded non-negativity 
Model = [replace(p,x,x0)==y0,
         replace(p,x,x1)>=y1, 
         replace(p,x,x2)==y2,
         replace(dp,x,x0)==0,
         replace(dp,x,x2)==0,
         replace(dp2,x,x1)>=0];
xgrid = linspace(-1,1,15);     
for i = 1:length(xgrid)
 Model = [Model, replace(p,x,xgrid(i)) >= 0];
end
optimize(Model, int(p^2,x,-1,1));
yv = polyval(fliplr(value(a')),xv);
ygrid = polyval(fliplr(value(a')),xgrid);
hold on
plot(xv,yv);
plot(xgrid,ygrid,'r+');
grid on

%% Non-negativity via SOS
Model = [replace(p,x,x0)==y0,
         replace(p,x,x1)>=y1, 
         replace(p,x,x2)==y2,
         replace(dp,x,x0)==0,
         replace(dp,x,x2)==0,
         replace(dp2,x,x1)>=0];
[s,c] = polynomial(x,2);
Model = [Model, sos(s), sos(p - s*(1-x^2))];
solvesos(Model, int(p^2,x,-1,1),[],[a;c]);
yv = polyval(fliplr(value(a')),xv);
hold on
plot(xv,yv,'--b');

%% Less conservative SOS model
Model = [replace(p,x,x0)==y0,
         replace(p,x,x1)>=y1,
         replace(p,x,x2)==y2,
         replace(dp,x,x0)==0,
         replace(dp,x,x2)==0,
         replace(dp2,x,x1)>=0];
[s,c] = polynomial(x,6);
Model = [Model, sos(s), sos(p - s*(1-x^2))];
solvesos(Model, int(p^2,x,-1,1),[],[a;c]);
yv = polyval(fliplr(value(a')),xv);
hold on
plot(xv,yv,'-b');