% % sparse_int.m % clear all % % Input constants % % % Compute the function on a uniform grid. % N1 = 41; N2 = 41; xvals1 = linspace(0,1,N1); xvals2 = linspace(0,1,N2); unit1 = ones(size(xvals1)); unit2 = ones(size(xvals2)); unit = ones(N1,N2); [X1,X2] = meshgrid(xvals1,xvals2); Z = X2.*exp(X1.*X2); % % Compute the interpolant % ell = 3; for j = 1:ell if j==1; n_ell = 1; qvals = 0; else n_ell = 2^(j-1)+1; N = n_ell-1; qvals = clencurt(N); end q = (1/2)*(1-qvals); if j==1 n1 = n_ell; x1 = q; y1 = q'; f0 = fun_exp(x1,y1,n1,n1); F01 = f0*ones(size(xvals1)); F02 = f0*ones(size(xvals2)); F0 = f0*ones(N1,N2); Iprev = f0*ones(N1,N2); elseif j==2 x2 = q; y2 = q'; n2 = n_ell; fvals12 = fun_exp(x1,y2,n1,n2); fvals21 = fun_exp(x2,y1,n2,n1); den1 = (x2(1)-x2(2))*(x2(1)-x2(3)); den2 = (x2(2)-x2(1))*(x2(2)-x2(3)); den3 = (x2(3)-x2(1))*(x2(3)-x2(2)); I12 = fvals12(1)*(xvals2-x2(2)).*(xvals2-x2(3))/den1 ... + fvals12(2)*(xvals2-x2(1)).*(xvals2-x2(3))/den2 ... + fvals12(3)*(xvals2-x2(1)).*(xvals2-x2(2))/den3 - F02; I21 = fvals21(1)*(xvals1-x2(2)).*(xvals1-x2(3))/den1 ... + fvals21(2)*(xvals1-x2(1)).*(xvals1-x2(3))/den2 ... + fvals21(3)*(xvals1-x2(1)).*(xvals1-x2(2))/den3 - F01; I = Iprev + unit1'*I12 + I21'*unit2; Iprev = I; else x3 = q; y3 = q'; n3 = n_ell; fvals13 = fun_exp(x1,y3,n1,n3); fvals31 = fun_exp(x3,y1,n3,n1); fvals22 = fun_exp(x2,y2,n2,n2); f12 = fvals22(1,2); f32 = fvals22(3,2); f21 = fvals22(2,1); f23 = fvals22(2,3); f22 = fvals22(2,2); den1 = (x3(1)-x3(2))*(x3(1)-x3(3))*(x3(1)-x3(4))*(x3(1)-x3(5)); den2 = (x3(2)-x3(1))*(x3(2)-x3(3))*(x3(2)-x3(4))*(x3(2)-x3(5)); den3 = (x3(3)-x3(1))*(x3(3)-x3(2))*(x3(3)-x3(4))*(x3(3)-x3(5)); den4 = (x3(4)-x3(1))*(x3(4)-x3(2))*(x3(4)-x3(3))*(x3(4)-x3(5)); den5 = (x3(5)-x3(1))*(x3(5)-x3(2))*(x3(5)-x3(3))*(x3(5)-x3(4)); I13 = fvals13(1)*(xvals2-x3(2)).*(xvals2-x3(3)).*(xvals2-x3(4)).*(xvals2-x3(5))/den1 ... + fvals13(2)*(xvals2-x3(1)).*(xvals2-x3(3)).*(xvals2-x3(4)).*(xvals2-x3(5))/den2 ... + fvals13(3)*(xvals2-x3(1)).*(xvals2-x3(2)).*(xvals2-x3(4)).*(xvals2-x3(5))/den3 ... + fvals13(4)*(xvals2-x3(1)).*(xvals2-x3(2)).*(xvals2-x3(3)).*(xvals2-x3(5))/den4 ... + fvals13(5)*(xvals2-x3(1)).*(xvals2-x3(2)).*(xvals2-x3(3)).*(xvals2-x3(4))/den5; I31 = fvals31(1)*(xvals1-x3(2)).*(xvals1-x3(3)).*(xvals1-x3(4)).*(xvals1-x3(5))/den1 ... + fvals31(2)*(xvals1-x3(1)).*(xvals1-x3(3)).*(xvals1-x3(4)).*(xvals1-x3(5))/den2 ... + fvals31(3)*(xvals1-x3(1)).*(xvals1-x3(2)).*(xvals1-x3(4)).*(xvals1-x3(5))/den3 ... + fvals31(4)*(xvals1-x3(1)).*(xvals1-x3(2)).*(xvals1-x3(3)).*(xvals1-x3(5))/den4 ... + fvals31(5)*(xvals1-x3(1)).*(xvals1-x3(2)).*(xvals1-x3(3)).*(xvals1-x3(4))/den5; den1 = (x2(1)-x2(2))*(x2(1)-x2(3)); den2 = (x2(2)-x2(1))*(x2(2)-x2(3)); den3 = (x2(3)-x2(1))*(x2(3)-x2(2)); I22 = fvals22(1,1)*((xvals1-x2(2)).*(xvals1-x2(3)))'*((xvals2-x2(2)).*(xvals2-x2(3)))/(den1*den1) ... + fvals22(2,1)*((xvals1-x2(1)).*(xvals1-x2(3)))'*((xvals2-x2(2)).*(xvals2-x2(3)))/(den2*den1) ... + fvals22(3,1)*((xvals1-x2(1)).*(xvals1-x2(2)))'*((xvals2-x2(2)).*(xvals2-x2(3)))/(den3*den1) ... + fvals22(1,2)*((xvals1-x2(2)).*(xvals1-x2(3)))'*((xvals2-x2(1)).*(xvals2-x2(3)))/(den1*den2) ... + fvals22(2,2)*((xvals1-x2(1)).*(xvals1-x2(3)))'*((xvals2-x2(1)).*(xvals2-x2(3)))/(den2*den2) ... + fvals22(3,2)*((xvals1-x2(1)).*(xvals1-x2(2)))'*((xvals2-x2(1)).*(xvals2-x2(3)))/(den3*den2) ... + fvals22(1,3)*((xvals1-x2(2)).*(xvals1-x2(3)))'*((xvals2-x2(1)).*(xvals2-x2(2)))/(den1*den3) ... + fvals22(2,3)*((xvals1-x2(1)).*(xvals1-x2(3)))'*((xvals2-x2(1)).*(xvals2-x2(2)))/(den2*den3) ... + fvals22(3,3)*((xvals1-x2(1)).*(xvals1-x2(2)))'*((xvals2-x2(1)).*(xvals2-x2(2)))/(den3*den3) ... - unit1'*(I12+F02) - (I21+F01)'*unit2; - 2*f22*unit - f12*unit - f32*unit - f21*unit - f23*unit; %I = Iprev + I31'*unit2 + unit1'*I13 + I22; I = I31'*unit2 + unit1'*I13 + I22; end end %I = Iprev; xvals = 1/2*ones(1,5); nodes1 = [x3' ; xvals]; nodes2 = nodes1'; nodes3 = [0 0.5 1 0 0.5 1 0 0.5 1; 0 0 0 0.5 0.5 0.5 1 1 1]; nodes = [0 0.5 1 0 0 0 0.5 0 x3(2) 0.5 x3(4) 1 0.5 0 0.5 1; 0 0 0 0 0.5 1.0 x3(2) 0.5 0.5 0.5 0.5 0.5 x3(4) 1 1 1]; xnodes = nodes(1,:); ynodes = nodes(2,:); znodes = ynodes.*exp(xnodes.*ynodes); % % Plot the function and observations % figure(1) contour(X1,X2,Z,20,'linewidth',2); axis('equal') hold on scatter(nodes1(1,:),nodes1(2,:),'b','linewidth',4) scatter(nodes1(2,:),nodes1(1,:),'b','linewidth',4) scatter(nodes3(1,:),nodes3(2,:),'b','linewidth',4) hold off set(gca,'Fontsize',[22]); xlabel('q_1') ylabel('q_2') %print -depsc full_contour figure(2) mesh(X1,X2,Z) set(gca,'Fontsize',[22]); xlabel('q_1') ylabel('q_2') figure(3) contour(X1',X2',I,20,'linewidth',2); axis('equal') hold on scatter(nodes1(1,:),nodes1(2,:),'b','linewidth',4) scatter(nodes1(2,:),nodes1(1,:),'b','linewidth',4) scatter(nodes3(1,:),nodes3(2,:),'b','linewidth',4) hold off set(gca,'Fontsize',[22]); xlabel('q_1') ylabel('q_2') %print -depsc sparse_contour figure(4) mesh(X1',X2',I) axis([0 1 0 1 0 3]); hold on plot3(xnodes,ynodes,znodes,'ok','linewidth',3) hold off set(gca,'Fontsize',[22]); xlabel('q_1') ylabel('q_2') zlabel('f(q_1,q_2)') %print -depsc function_sparse_pts int_points = [I(1,1) I(1,21) I(1,41) I(21,1), I(21,21), I(21,41), I(41,1), I(41,21), I(41,41)] Z_points = [Z(1,1) Z(1,21) Z(1,41) Z(21,1), Z(21,21), Z(21,41), Z(41,1), Z(41,21), Z(41,41)]