Can contourplot plot based on the peaks of the function?

Hi,

I need to plot a contourplot based on the peaks of my function. I know the maximum of my function but the contourplot maximum is different than that. I need to plot something like the below fig. (the data in 2 figures are not the same)

but I get this one in which, the real peak and the contour peak are not the same. Is there any command for edditing this?

Thank.

[Cris LaPierre]

Share your code, or at a minimum the data.

Answer 1

[Cris LaPierre]

I suspect there is an error in your code. Here's a simple example using the peaks funciton in MATLAB

[x,y,z] = peaks(100);
surf(x,y,z)
view(2)

figure
contour(x,y,z,20)

Here is what the peak function would look like from the top down (view(2))

And here is what the corresponding contour plot looks like. Notice that the peaks align.

[Cris LaPierre]

Since it may be a little tough to visualize in 2D, here's a 3D plot of the same data. Three peaks that go up, 2 peaks that go down.

[ASH]

Thanks for your reply,

My variables and function are more complex than simple vectors and I'm not sure how to use peaks for them. Moreover, I think because one of my variables is a large number I got error of "maximum array size" when using meshgrid. Also I need to have three of these contour plots (for another code) in one figure, but "hold on" doesn't work.

Because it's a long code I preferred not to share it but if you need it I attached that below. Hopefully, you can find out what's going on.

clear
clc

global T
global K
global pbs
global V_bpi


T=[1600;1600;1600;1600;1600;2000;2000;2000;2000;2000];
K=[2.72925e9;1.43158e9;2.41354e9;2.17249e9;3.89502e8;1.05276e10;1.03348e10;1.08849e10;1.07148e10;1.07487e10];
pbs=[6.19204e12;14354];
psd=[6.03563e12;514.848];
n=length(T);
p=2;

%best estimates
fun = @(x)sseval(x,T,K);
x0 = [1.2e13,14000];
bestx = fminsearch(fun,x0);
c=[bestx(1) bestx(2)]; 

K_p=c(1)*exp(-c(2)./T);
err=sum((K_p-K).^2)/(n-p);

f=exp(-c(2)./T);
g=-(c(1)./T).*exp(-c(2)./T);
F=[f g];
Ft=F';
FtF=Ft*F;

%covariance matrix
V_br=inv(FtF)*err;
cor=V_br(1,2)/(V_br(1,1)*V_br(2,2))^.5;
V_bp=[psd(1),cor*(psd(1)*psd(2))^.5;cor*(psd(1)*psd(2))^.5,psd(2)];
V_bpi=inv(V_bp);
det=det(V_bp);


%Correlated
P_ab=@Calculate_P_ab;
Mat=@Matrix;

D=@(te1,te2)exp(-.5*P_ab(te1,te2)).*exp(-Mat(te1,te2)/(2*err))/((2*pi*err)^(n/2)*((2*pi)^2*det)^.5);
% Y_p=integral2(D,1.1e13,1.2e13,14000,14100);
P_y=@(te1,te2)D(te1,te2);
x0=[7e12 13900];
[xmax,fval]=fmincon(@(x)-P_y(x(1),x(2)),x0);

% figure
% surf(J,H,D(J,H), 'EdgeColor','none', 'FaceAlpha',0.5)
[J,H] = meshgrid(1.1e13:1.2e13,14000:14200);

figure
contour(J,H,P_y(J,H),10)
hold on

plot3(xmax(1),xmax(2), P_y(xmax(1),xmax(2)), 'gp', 'MarkerSize', 10, 'MarkerFaceColor','g')
hold off
xlabel('slope')
ylabel('intercept')

function val_bs = sseval(x,T,K)
  A=x(1);
  ER=x(2);
    val_bs=sum(K-(A*exp(-ER./T))).^2;  
end

function Val=Matrix(te1,te2)
    Val=0;
    global T
    global K
    global n
    for k=1:n 
        var1=T(k);
        var2=K(k);
        Val=Val+(var2-(te1.*exp(-te2./var1))).^2;
    end
end


function val_ab = Calculate_P_ab(te1,te2)
    global V_bpi
    global pbs
    val_ab=V_bpi(1,1)*(te1-pbs(1)).^2+2*V_bpi(1,2)*(te1-pbs(1)).*(te2-pbs(2))+V_bpi(2,2)*(te2-pbs(2)).^2;
end

Thanks.

[Cris LaPierre]

Don't use peaks. You hadn't shared your data, so I used peaks to create some data so I could then generate a contour plot.

[Cris LaPierre]

When I tried to run your code, all the values of P_y(J,H) were zero. Perhaps start by making sure your calculations are correct? If there is no variation in Z, you don't want to use contour.

I don't know how many points you need in your meshgrid, but 1e12 seems like too many. Start small and increase only if necessary. I tried this before encountering the issue with P_y.

[J,H] = meshgrid(linspace(1.1e13,1.2e13,1000),linspace(14000,14200));

[ASH]

Thanks again for your response,

I corrected my code and I got good results like the below fig. I wonder how you read the final value of a function handle (P_y here) because I need that for another code. I tried "pretty" and it doesn't work.

Thanks

[Walter Roberson]

If you have the symbolic toolbox, then

syms Te1 Te2
P_y(Te1, Te2)

and you might also want to simplify() that.

[Cris LaPierre]

I got the values by doing this:

[J,H] = meshgrid(linspace(1.1e13,1.2e13,1000),linspace(14000,14200));
a=P_y(J,H);

I then opened a in the variable editor (double click on it in the workspace).

[ASH]

Thanks a lot!