Why the following code is not generating plot for W?

1 次查看(过去 30 天)
Reference figure: plot should be like this for W refer square ones.
Ho = 1;
alpha = 0.1;
a = 1.0;
sigma = 0.15;
lbar_list = [0.0, 0.2];
eps = 0:0.2;
hold on
for i = 1:numel(lbar_list)
lbar = lbar_list(i);
H = @(x) Ho + a*(1 - x);
G1 = @(x,eps) H(x).^3 + 3 .* H(x).^2 .* alpha + 3 .* H(x) .* alpha^2 + 3 .* H(x) .* sigma^2 + eps + 3*sigma^2*alpha + alpha^3 - 12.*lbar^2 .* (H(x) + alpha);
G2 = @(x) 24 .* lbar.^3 .* tanh(H(x)./(2.*lbar));
G3 = @(x,eps) (12.*lbar^2*alpha - eps - alpha^3 - 3*sigma^2*alpha) .* (1 - (tanh(H(x)./(2.*lbar))).^2);
G = @(x,eps) G1(x,eps) + G2(x) + G3(x,eps);
Hm1 = @(x,eps) H(x).* (1 ./ G(x,eps));
Hm2 = @(x,eps) (1 ./ G(x,eps));
IntHm1 =@(eps) integral(@(x) Hm1(x,eps),0,1);
IntHm2 =@(eps) integral(@(x) Hm2(x,eps),0,1);
Hm = @(eps) IntHm1(eps) / IntHm2(eps);
P1 = @(x,eps) 6 .* (1 ./ G(x,eps)) .* (H(x) - Hm(eps));
P2 = @(x,eps) integral(@(x) P1(x,eps),0,x);
W = @(eps) integral(@(x) P2(x,eps),0,1, 'ArrayValued', true);
plot(eps,4.5*W(eps))
end
legend(num2str(lbar_list'))
ylim([0.4 0.8])
xlim([0 0.2])
set(gca, 'ytick', 0.4:0.1:0.8);
set(gca, 'xtick', 0:0.02:0.2);
xlabel('eps')
ylabel('W')

采纳的回答

KSSV
KSSV 2022-7-16
I remember you asking the similiar question earlier. What have you learned from the previous question? You are keep on asking same questions. Spend some time on the documentation and rethink on your old codes.
clc; clear all ;
Ho = 1;
alpha = 0.1;
a = 1.0;
sigma = 0.15;
lbar_list = [0.0, 0.2];
eps = linspace(0,0.2);
hold on
for i = 1:numel(lbar_list)
lbar = lbar_list(i);
H = @(x) Ho + a*(1 - x);
G1 = @(x,eps) H(x).^3 + 3 .* H(x).^2 .* alpha + 3 .* H(x) .* alpha^2 + 3 .* H(x) .* sigma^2 + eps + 3*sigma^2*alpha + alpha^3 - 12.*lbar^2 .* (H(x) + alpha);
G2 = @(x) 24 .* lbar.^3 .* tanh(H(x)./(2.*lbar));
G3 = @(x,eps) (12.*lbar^2*alpha - eps - alpha^3 - 3*sigma^2*alpha) .* (1 - (tanh(H(x)./(2.*lbar))).^2);
G = @(x,eps) G1(x,eps) + G2(x) + G3(x,eps);
Hm1 = @(x,eps) H(x).* (1 ./ G(x,eps));
Hm2 = @(x,eps) (1 ./ G(x,eps));
IntHm1 =@(eps) integral(@(x) Hm1(x,eps),0,1);
IntHm2 =@(eps) integral(@(x) Hm2(x,eps),0,1);
Hm = @(eps) IntHm1(eps) / IntHm2(eps);
P1 = @(x,eps) 6 .* (1 ./ G(x,eps)) .* (H(x) - Hm(eps));
P2 = @(x,eps) integral(@(x) P1(x,eps),0,x);
W = @(eps) integral(@(x) P2(x,eps),0,1, 'ArrayValued', true);
WW = zeros(1,length(eps)) ;
for j = 1:length(eps)
WW(j) = 4.5*W(eps(j)) ;
end
plot(eps,WW)
end
legend(num2str(lbar_list'))
ylim([0.4 0.8])
xlim([0 0.2])
set(gca, 'ytick', 0.4:0.1:0.8);
set(gca, 'xtick', 0:0.02:0.2);
xlabel('eps')
ylabel('W')

更多回答(0 个)

类别

Help CenterFile Exchange 中查找有关 Particle & Nuclear Physics 的更多信息

Community Treasure Hunt

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

Start Hunting!

Translated by