Clear Filters
Clear Filters

Why the following code is not generating plot for W?

1 view (last 30 days)
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')

Accepted Answer

KSSV
KSSV on 16 Jul 2022
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')
  1 Comment
Torsten
Torsten on 16 Jul 2022
@AVINASH SAHU comment moved here:
Okay thanks I will spend more time reading the documentation

Sign in to comment.

More Answers (0)

Categories

Find more on Particle & Nuclear Physics in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!