approx_int_func_solve.m

function ww = approx_int_func_solve(xy_pairs, xmin, xmax, intval, w0, f, lbda)
% example: 
% xy_pairs = [0, 1.32; 33.9, .55];
% xmin = 0;
% xmax = 33.9;
% intval = 38.985;
% w0 = [0.1,-2,2]
% ww = approx_int_func_solve(xy_pairs, xmin, xmax, intval, w0);

if nargin < 6
  f = @(x,w) -exp(w(1)*x+w(2)) + w(3);
  lbda = 0.01;
end

% \sum_i [y_i - f(x_i;w)]^2 + lbda * [\int_xmin^xmax f(x;w) dx - intval]^2
L = @(w) norm(xy_pairs(:,2) - arrayfun(@(x) f(x,w), xy_pairs(:,1)), 2).^2 + lbda * (integral(@(x) f(x,w), xmin, xmax) - intval)^2;

options = optimoptions(@fminunc, 'Algorithm', 'quasi-newton', 'Display', 'iter');
[ww,Lval,exitflag,output] = fminunc(L, w0, options);