Please find the mistake in this code.
2 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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 commentaires
Sam Chak
le 8 Juil 2022
Can you plot the lines manually one by one? You'll definitely find out what went wrong then.
Réponse acceptée
Torsten
le 8 Juil 2022
Modifié(e) : Torsten
le 8 Juil 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)')
Plus de réponses (0)
Voir également
Catégories
En savoir plus sur Logical dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!

