# Why the following code is not generating plot for W?

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

### Categories

Find more on Histograms 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!