How to plot a 3D surface from scatter data inside a loop?

1 次查看(过去 30 天)
Hello there, this is my first post here, and i'd like to thank you all in advance. So my problem is that i'm trying to plot a 3D surface from the scatter data result of an ODE inside a loop. I'm using this to understand how this variables interfere in a electric field used in a quantum control scenario. The thing is that the relation between the variables is not direct, ie, the X and Y components are inside a loop and they are used to create an expression that is interpolated and than used to solve the ODE. I can plot the individual points using Scatter3, but i can't use those points to make a surface. The code i'm using is below. How should i proceed? Thanks again!
% Variables
ft = -20000:.1:20000;
for alpha = 0.001:0.0005:0.01
for beta = 0.001:0.0005:0.01
ai = 0.4;
ai2 = 0.6;
af = 1;
w0 = 0.02;
mi = 6;
phii = 0;
phif = pi/3.2;
% Functions
g = 1./(1+exp(-alpha.*ft));
f = ai*(1-g)+af*g;
p = 1./(1+exp(-beta.*ft));
h = phii*(1-p)+phif*p;
% Electric Field
E = alpha.*(af-ai).*exp(alpha.*ft).*sin(w0.*ft+h)./(mi.*(1+exp(alpha.*ft)).*sqrt((1-ai+(1-af).*exp(alpha.*ft)).*(ai+af.*exp(alpha.*ft))))+2.*beta.*(phif-phii).*exp(beta.*ft).*sqrt(f.*(1-f)).*cos(w0.*ft+h)/(mi.*(1-2.*f).*(1+exp(beta.*ft).^2));
% ODE
[t,y] = ode45(@(t,y) myode(t,y,ft,E,mi,w0), ft, [sqrt(ai)*exp(1i.*0) sqrt(ai2)*exp(1i.*0)]);
z=1-abs(y(400000,1)).^2;
% Plot
scatter3(alpha,beta,z)
xlabel('Alpha')
ylabel('Beta')
hold on
N = 50;
xi = linspace(min(alpha),max(alpha),N);
yi = linspace(min(beta),max(beta),N);
[X,Y] = meshgrid(xi,yi);
Z = griddata(alpha,beta,z,X,Y);
surf(X,Y,Z)
end
end
% End of the program
function dydt = myode(t,y,ft,E,mi,w0)
E = interpn(ft,E,t);
dydt = [1i.*mi.*y(2).*E.*exp(-1i.*w0.*t);1i.*mi.*y(1).*E.*exp(1i.*w0.*t)];
end

采纳的回答

Star Strider
Star Strider 2019-4-8
Your data are gridded, however you need to use the reshape function to create matrices from them. I changed your code slightly to save the relevant variables to a matrix (after preallocating ‘abz’ before the loop, and introducing a counter variable ‘kntr’ at the start of the ‘beta’ loop, and took the scatter3 call completely out of the loops):
for beta = 0.001:0.0005:0.01
kntr = kntr + 1
then:
% ODE
[t,y] = ode45(@(t,y) myode(t,y,ft,E,mi,w0), ft, [sqrt(ai)*exp(1i.*0) sqrt(ai2)*exp(1i.*0)]);
z=1-abs(y(400000,1)).^2;
abz(kntr,:) = [alpha, beta, z];
and then after the loops:
ABZ = table(abz(:,1), abz(:,2), abz(:,3), 'VariableNames',{'alpha','beta','z'})
writetable(ABZ, fullfile(dirpath, 'ABZ_20190408.txt'))
with ‘dirpath’ being the directory you want to write the file to. I attached the file to my Answer.
The plot is then created as:
ABZ = readtable('ABZ_20190408.txt');
alphamtx = reshape(ABZ.alpha, 19, []);
betamtx = reshape(ABZ.beta, 19, []);
zmtx = reshape(ABZ.z, 19, []);
N = 50;
xi = linspace(min(ABZ.alpha),max(ABZ.alpha),N);
yi = linspace(min(ABZ.beta),max(ABZ.beta),N);
[X,Y] = meshgrid(xi,yi);
Z = griddata(alphamtx,betamtx,zmtx,X,Y);
figure
surfc(X,Y,Z)
grid on
xlabel('\alpha')
ylabel('\beta')
zlabel('z')
shading('interp')
And your complete revised code (except for the table and writetable calls) is:
% Variables
kntr = 0;
abz = zeros(19^2, 3);
ft = -20000:.1:20000;
for alpha = 0.001:0.0005:0.01
for beta = 0.001:0.0005:0.01
kntr = kntr + 1
ai = 0.4;
ai2 = 0.6;
af = 1;
w0 = 0.02;
mi = 6;
phii = 0;
phif = pi/3.2;
% Functions
g = 1./(1+exp(-alpha.*ft));
f = ai*(1-g)+af*g;
p = 1./(1+exp(-beta.*ft));
h = phii*(1-p)+phif*p;
% Electric Field
E = alpha.*(af-ai).*exp(alpha.*ft).*sin(w0.*ft+h)./(mi.*(1+exp(alpha.*ft)).*sqrt((1-ai+(1-af).*exp(alpha.*ft)).*(ai+af.*exp(alpha.*ft))))+2.*beta.*(phif-phii).*exp(beta.*ft).*sqrt(f.*(1-f)).*cos(w0.*ft+h)/(mi.*(1-2.*f).*(1+exp(beta.*ft).^2));
% ODE
[t,y] = ode45(@(t,y) myode(t,y,ft,E,mi,w0), ft, [sqrt(ai)*exp(1i.*0) sqrt(ai2)*exp(1i.*0)]);
z=1-abs(y(400000,1)).^2;
abz(kntr,:) = [alpha, beta, z];
end
end

更多回答(0 个)

类别

Help CenterFile Exchange 中查找有关 Historical Contests 的更多信息

产品

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by