examp.m

这是在网上下的一个东东

%
% Reset to a known state.
%
clear
rand('state',0);
randn('state',0);
%
%
%
disp('This takes about 5 minutes on a 3 Ghz Pentium 4');
%
% This script creates the synthetic data set for the EM-38 example.
%
%
% Number of layers.
%
M=11;
%
% The model conductivities in S/m.
%
mtrue=[100; 95; 90; 95; 100; 130; 160; 200; 250; 310; 360]/1000.0;
%
% Layer thicknesses.
%
D=0.2*ones(10,1);
%
% Operating frequency.
%
f=14600;
OMEGA=2*pi*f;
%
% Distance between coils.
%
R=1.0;
%
% Conductivity of the air above the top layer.
%
SIGMA0=0;
%
% Magnetic permeabilities.
%
MU0=pi*4e-7;
MU=MU0*ones(11,1);
%
% Scaling factor Delta.
%
DELTA=sqrt(2/(mtrue(1)*MU0*OMEGA));
%
% Heights of measurements.
%
heights=[0.0; 0.1; 0.2; 0.3; 0.4; 0.5; 0.75; 1.0; 1.50];
%
% Now, do the predictions.
%
pred=[];
pred=zeros(length(heights),2);
for i=1:length(heights)
  H=heights(i);
  [predv,predh]=predict(R,DELTA,H,M,MU,MU0,mtrue,SIGMA0,D,OMEGA);
  pred(i,:)=[predv predh];
end
%
% Add random noise.
%
rand('state',0);
datanf=[pred(:,1); pred(:,2)];
data=datanf+0.1*randn(size(datanf));
%
% Save the data set.
%
%save data.mat data 
%
% Now, solve the problem.
%
m0=200*ones(11,1)/1000;
%
% First, try using LM.
%
global DATA;
DATA=data;
mlm=lm('fund','jac',m0,1.0e-6,50);
disp('The chi^2 value is ');
norm(fund(mlm),2)^2
disp('With 18-11=9 dof');
disp('cond(J''*J) is ');
Jlm=jac(mlm);
cond(Jlm'*Jlm)
%
% Plot the LM solution.
%
figure(1);
bookfonts
plotconst(mlm*1000,0,2.2);
xlabel('depth (m)');
ylabel('Conductivity (mS/m)');
%print -deps mlm.eps
%
% Next, occam.
%
L2=full(get_l(11,2));
m0=200*ones(11,1)/1000;
delta=0.1*sqrt(18);
moccam=occam('fun','jac',L2,data,m0,delta);
figure(2);
bookfonts;
plotconst(moccam*1000,0,2.2);
xlabel('depth (m)');
ylabel('Conductivity (mS/m)');
%print -deps moccam.eps
%
% Plot the true model.
%
figure(3)
bookfonts;
plotconst(mtrue*1000,0,2.2);
xlabel('depth (m)');
ylabel('Conductivity (mS/m)');


%print -deps mtrue.eps
%
% Save the whole solution.
%
%save examp.mat