u/johannjc/injection: test_solver2.m

File test_solver2.m, 1.4 KB (added by Jonathan, 7 years ago)
Line 
1% initialize data with ideal solution
2x0=.05
3y0=.05
4[x,y]=meshgrid(linspace(-1,1,10), linspace(-1,1,10))
5dv=(x(1,2)-x(1,1)).^2;
6mass_target=pi*(1-exp(-1))
7rho_predict=arrayfun(@(x,y) exp(-((x-x0)^2+(y-y0)^2)).*((x-x0).^2+(y-y0).^2<1), x, y)/mass_target;
8mask=((x-x0).^2+(y-y0).^2 < 1);
9
10
11% Plot initial guess
12subplot(1,4,1)
13imagesc(x(1,:),y(:,1),rho_predict)
14colorbar()
15
16%Predict correct mass, and center
17mass_predict=sum(sum(rho_predict))*dv
18com_x=sum(sum(rho_predict.*x))*dv/mass_predict
19com_y=sum(sum(rho_predict.*y))*dv/mass_predict
20title(['initial guess - com = (',num2str(com_x),',',num2str(com_y),')', ' mass=', num2str(mass_predict)])
21
22
23% Calculate solution in least squares sense
24A=[x(mask)';y(mask)';ones(1,numel(x(mask)))]*dv;
25B=[x0; y0; mass_target];
26sol=pinv(A)*B;
27rho_correct=zeros(size(rho_predict));
28rho_correct(mask)=sol;
29
30subplot(1,4,2)
31imagesc(x(1,:), y(:,1), rho_correct)
32colorbar()
33title('Solution with smallest norm')
34
35sol=lsqlin(eye(numel(x(mask))), rho_predict(mask), [],[],A,B)
36subplot(1,4,3)
37rho_correct=zeros(size(rho_predict));
38rho_correct(mask)=sol;
39imagesc(x(1,:), y(:,1),rho_correct)
40title('Solution closest to prediction')
41colorbar()
42
43
44sol=rho_predict(mask)-A'*inv(A*A')*(A*rho_predict(mask) - B);
45subplot(1,4,4)
46rho_correct=zeros(size(rho_predict));
47rho_correct(mask)=sol;
48imagesc(x(1,:), y(:,1),rho_correct)
49title('Solution closest to prediction')
50colorbar()