examp.m

这是在网上下的一个东东

%
% This script generates the G matrix, mtrue, dtrue, and noisy d for the
% source history example.   
%
clear 
rand('state',0);
randn('state',0);
% set up input for domre.m
params=[1.0, 1.0];   % [v,D]
xsample=[0.01,25.05,50:10:250,275,300];
tsample=300;
nt=100;
tmax=250.;
tmin=0.01;
t=linspace(tmin,tmax,nt);
%
% get true source concentration as a function of t
%
cin=cintrue(t);
%
% Make the G matrix.
%
G=getrep('kernel',t,xsample,tsample,params);
% generate sampled data
cexact=G*cin';
rand('state',0);
csample=cexact+0.001*randn(size(cexact));
%
%
%
mtrue=cin';
d=csample;
dtrue=cexact;
%save sourcehist.mat mtrue d G dtrue xsample t
%
% Plot out the cin and data.
%
figure(1);
bookfonts
plot(t,cin,'ok-');
xlabel('Time');
ylabel('C_{in}(t)');
%print -deps cintrue.eps
%
% Plot the downstream samples.
%
figure(2);
bookfonts
plot(xsample,d,'ok-');
xlabel('Distance');
ylabel('C(x,300)');
%print -deps cout.eps
%
% Save the problem.
%
%save sourcehist.mat
%
% This script uses several different methods to solve the source history
% reconstruction problem.
%
%
% Get the problem size.
%
[m,n]=size(G);
%
% First, get the least squares solution.
%%
mls=pinv(G)*d;
figure(3);
bookfonts
plot(t,mls,'ok-');
ylabel('C_{in}(t)');
xlabel('Time');
%print -deps mls
%
% Next, solve the problem using lsqnonneg
%
mlsn=lsqnonneg(G,d);
figure(4);
bookfonts
plot(t,mlsn,'ok-');
xlabel('Time');
ylabel('C_{in}(t)');
%print -deps mlsn
%
% Next, solve the problem using BVLS
%
l=zeros(n,1);
u=1.1*ones(n,1);
mbvls=bvls(G,d,l,u);
figure(5);
bookfonts
plot(t,mbvls,'ok-');
xlabel('Time');
ylabel('C_{in}(t)');
%print -deps mbvls
%
% Next, solve the problem using tikhcstr
%
lambdas=logspace(-6,1,20);
rhos=zeros(size(lambdas));
eta=zeros(size(lambdas));
L=full(get_l(n,2));
A=[eye(n); -eye(n)];
b=[l; -u];
for i=1:length(lambdas)
  [m,rho(i),eta(i)]=tikhcstr(G,d,A,b,L,lambdas(i));
end  
figure(6);
bookfonts
loglog(rho,eta,'ok-');
xlabel('|| Gm-d ||');
ylabel('|| Lm ||');
%
% Put in dotted lines showing the L-curve corner.
%
hold on
plot([0.001 rho(14)],[eta(14) eta(14)],'k:');
plot([rho(14) rho(14)],[0.01 eta(14)],'k:');
%print -deps lcurve.eps
%
% Get the corner model and plot it.
%
lambda=lambdas(14);
mtik=tikhcstr(G,d,A,b,L,lambda);
figure(7);
bookfonts
plot(t,mtik,'ok-');
xlabel('Time');
ylabel('Concentration');
%print -deps mtik
%
% Now, construct bounds on the minimum/maximum average value in the 
% interval from t=125 to t=150.  
%
c=zeros(n,1);
c(51:60)=0.1*ones(10,1);
[xmin,xmax]=blf2(G,d,c,0.01,l,u);
cminus=c'*xmin
cplus=c'*xmax
%
% Plot the xmin and xmax solutions.
%

figure(8);
bookfonts
plot(t,xmin,'ok-');
xlabel('Time');
ylabel('Concentration');
%print -deps xmin

figure(9);
bookfonts
plot(t,xmax,'ok-');
xlabel('Time');
ylabel('Concentration');
%print -deps xmax