Please find the mistake in this code.
3 views (last 30 days)
Show older comments
Reference image: plot should be like this but don't know where I am making mistake.
syms x lambda
warning off
alpha = -0.1;
sigma = 0.1;
eps = -0.1;
e = 0.2;
M = 10;
a = 2;
figure
psi_list = [10, 20, 30, 40. 50];
for i = 1:numel(psi_list)
psi = psi_list(i);
hbar = @(x) a - a.*x + x;
A1 = eps + alpha^3 + (3 * sigma^2 * alpha);
B1 =@(x,lambda) (-3 * lambda * M) * ((hbar(x).^2) + (2 .* hbar(x) .* alpha) + (sigma^2) + (alpha^2));
a1 =@(x) tanh(M .* hbar(x));
b1 =@(x) 1 - ((tanh(M .* hbar(x))).^2);
c1 = (M * alpha) - ((M^3 * A1)/3);
d1 = @(lambda) 2 * (M^2) * (1 + lambda);
C1=@(x) a1(x) + (b1(x) .* c1);
D1 =@(x,lambda) d1(lambda) .* ((hbar(x).^3) + (3 .* (hbar(x).^2) .* alpha) + (3 .* hbar(x) .* (alpha)^2) + (3 .* hbar(x) .* (sigma)^2) + eps + (3 * alpha * (sigma^2)) + (alpha^3));
f1 =@(x,lambda) B1(x,lambda) + (3 * lambda .* C1(x) .* hbar(x)) + (3 * lambda .* C1(x) .* alpha) + (D1(x,lambda) .* C1(x));
f2 =@(x,lambda) 12 * (M^2) * (1 + lambda) .* C1(x);
f3 = psi * (e^3);
f4 = @(lambda) (1 + lambda) *180 * ((1 - e)^2);
f5 =@(lambda) 1/(2 + lambda);
F = @(x,lambda) ((f5(lambda) .* f1(x,lambda))./f2(x,lambda)) + (f3/f4(lambda));
q1 =@(x,lambda) hbar(x) ./ (2 .* F(x,lambda));
Q1 = int(q1,x,0,1); % integration with respect to "x" limit 0 to 1
q2 = @(x,lambda) 1./(F(x,lambda));
Q2 = int(q2,x,0,1); % integration with respect to "x" limit 0 to 1
Q = @(lambda) Q1(lambda)/Q2(lambda); % Q is a function of "lambda"
p1 = @(x,lambda) (1./F(x,lambda)) .* ((0.5 .* hbar(x)) - Q(lambda));
P = @(x,lambda) int(p1,x); % integration with respect to "x" limit 0 to x
% W is obtained by integrating "P" with respect to x limit to x is 0 to
% 1
W =@(lambda) int(P,x,0,1);
fplot(@(lambda) W(lambda))
end
legend(num2str(psi_list'))
ylim([0 1])
xlim([0 4])
set(gca, 'ytick', 0:0.1:1);
set(gca, 'xtick', 0:0.2:4);
xlabel('lambda')
ylabel('W(lambda)')
6 Comments
Sam Chak
on 8 Jul 2022
Can you plot the lines manually one by one? You'll definitely find out what went wrong then.
Accepted Answer
Torsten
on 8 Jul 2022
Edited: Torsten
on 8 Jul 2022
alpha = -0.1;
sigma = 0.1;
eps = -0.1;
e = 0.2;
M = 10;
a = 2;
%figure
psi_list = [10, 20, 30, 40, 50];
lambda = 0:0.2:4;
for i = 1:numel(psi_list)
psi = psi_list(i);
hbar = @(x) a - a.*x + x;
A1 = eps + alpha^3 + (3 * sigma^2 * alpha);
B1 =@(x,lambda) (-3 * lambda * M) * ((hbar(x).^2) + (2 .* hbar(x) .* alpha) + (sigma^2) + (alpha^2));
a1 =@(x) tanh(M .* hbar(x));
b1 =@(x) 1 - ((tanh(M .* hbar(x))).^2);
c1 = (M * alpha) - ((M^3 * A1)/3);
d1 = @(lambda) 2 * (M^2) * (1 + lambda);
C1=@(x) a1(x) + (b1(x) .* c1);
D1 =@(x,lambda) d1(lambda) .* ((hbar(x).^3) + (3 .* (hbar(x).^2) .* alpha) + (3 .* hbar(x) .* (alpha)^2) + (3 .* hbar(x) .* (sigma)^2) + eps + (3 * alpha * (sigma^2)) + (alpha^3));
f1 =@(x,lambda) B1(x,lambda) + (3 * lambda .* C1(x) .* hbar(x)) + (3 * lambda .* C1(x) .* alpha) + (D1(x,lambda) .* C1(x));
f2 =@(x,lambda) 12 * (M^2) * (1 + lambda) .* C1(x);
f3 = psi * (e^3);
f4 = @(lambda) (1 + lambda) *180 * ((1 - e)^2);
f5 =@(lambda) 1./(2 + lambda);
F = @(x,lambda) ((f5(lambda) .* f1(x,lambda))./f2(x,lambda)) + (f3./f4(lambda));
q1 =@(x,lambda) hbar(x) ./ (2 .* F(x,lambda));
Q1 = @(lambda)integral(@(x)q1(x,lambda),0,1,'ArrayValued',true); % integration with respect to "x" limit 0 to 1
q2 = @(x,lambda) 1./(F(x,lambda));
Q2 = @(lambda)integral(@(x)q2(x,lambda),0,1,'ArrayValued',true); % integration with respect to "x" limit 0 to 1
Q = @(lambda) Q1(lambda)./Q2(lambda); % Q is a function of "lambda"
p1 = @(x,lambda) (1./F(x,lambda)) .* ((0.5 .* hbar(x)) - Q(lambda));
P = @(x,lambda) integral(@(y)p1(y,lambda),0,x,'ArrayValued',true); % integration with respect to "x" limit 0 to x
% W is obtained by integrating "P" with respect to x limit to x is 0 to
% 1
W =@(lambda) integral(@(x)P(x,lambda),0,1,'ArrayValued',true);
%fplot(@(lambda) W(lambda))
plot(lambda,W(lambda))
hold on
end
legend(num2str(psi_list'))
ylim([0 1])
xlim([0 4])
set(gca, 'ytick', 0:0.1:1);
set(gca, 'xtick', 0:0.2:4);
xlabel('lambda')
ylabel('W(lambda)')
More Answers (0)
See Also
Categories
Find more on Creating and Concatenating Matrices 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!