Why it is not making plot for different values of "M"?

1 次查看(过去 30 天)
syms x
alpha = -0.1;
sigma = 0.1;
eps = -0.1;
e = 0.2;
a = 2;
lambda = 2;
psi = 10;
figure
M_list = [4, 6, 8, 10];
for i = 1:numel(M_list)
M = M_list(i);
hbar = @(x) a - a.*x + x;
A1 = eps + alpha^3 + (3 * sigma^2 * alpha);
B1 = @(x) (-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 = 2 * (M^2) * (1 + lambda);
C1 = @(x) a1(x) + (b1(x) .* c1);
D1 = @(x) d1 .* ((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) B1(x) + (3 * lambda .* C1(x) .* hbar(x)) + (3 * lambda .* C1(x) .* alpha) + (D1(x) .* C1(x));
f2 = @(x) 12 * (M^2) * (1 + lambda) .* C1(x);
f3 = psi * (e^3);
f4 = (1 + lambda) *180 * ((1 - e)^2); % 180 is not given in paper
f5 = 1/(2 + lambda);
F = @(x) ((f5 .* f1(x))./f2(x)) + (f3/f4);
q1 = @(x) hbar(x) ./ (2 .* F(x));
Q1 = integral(q1,0,1);
q2 = @(x) 1./(F(x));
Q2 = integral(q2,0,1);
Q = Q1/Q2;
p1 = @(x) (1./F(x)) .* ((0.5 .* hbar(x)) - Q);
P = @(x) integral(p1,0,x);
fplot(P, [0 1])
ylim([0 1])
set(gca, 'ytick', 0:0.1:1);
set(gca, 'xtick', 0:0.2:1);
xlabel('x')
ylabel('P(x)')
end
Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is 1.6e+02. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.
Warning: Function behaves unexpectedly on array inputs. To improve performance, properly vectorize your function to return an output with the same size and shape as the input arguments.
Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is 8.0e+02. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.
Warning: Function behaves unexpectedly on array inputs. To improve performance, properly vectorize your function to return an output with the same size and shape as the input arguments.
Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is 1.7e+02. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.
Warning: Function behaves unexpectedly on array inputs. To improve performance, properly vectorize your function to return an output with the same size and shape as the input arguments.
Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is 2.1e+02. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.
Warning: Function behaves unexpectedly on array inputs. To improve performance, properly vectorize your function to return an output with the same size and shape as the input arguments.

采纳的回答

KSSV
KSSV 2022-7-8
编辑:KSSV 2022-7-8
You have to use hold on.
syms x
warning off
alpha = -0.1;
sigma = 0.1;
eps = -0.1;
e = 0.2;
a = 2;
lambda = 2;
psi = 10;
M_list = [4, 6, 8, 10];
figure
hold on
for i = 1:numel(M_list)
M = M_list(i);
hbar = @(x) a - a.*x + x;
A1 = eps + alpha^3 + (3 * sigma^2 * alpha);
B1 = @(x) (-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 = 2 * (M^2) * (1 + lambda);
C1 = @(x) a1(x) + (b1(x) .* c1);
D1 = @(x) d1 .* ((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) B1(x) + (3 * lambda .* C1(x) .* hbar(x)) + (3 * lambda .* C1(x) .* alpha) + (D1(x) .* C1(x));
f2 = @(x) 12 * (M^2) * (1 + lambda) .* C1(x);
f3 = psi * (e^3);
f4 = (1 + lambda) *180 * ((1 - e)^2); % 180 is not given in paper
f5 = 1/(2 + lambda);
F = @(x) ((f5 .* f1(x))./f2(x)) + (f3/f4);
q1 = @(x) hbar(x) ./ (2 .* F(x));
Q1 = integral(q1,0,1);
q2 = @(x) 1./(F(x));
Q2 = integral(q2,0,1);
Q = Q1/Q2;
p1 = @(x) (1./F(x)) .* ((0.5 .* hbar(x)) - Q);
P = @(x) integral(p1,0,x);
fplot(P, [0 1])
end
legend(num2str(M_list'))
ylim([0 1])
set(gca, 'ytick', 0:0.1:1);
set(gca, 'xtick', 0:0.2:1);
xlabel('x')
ylabel('P(x)')

更多回答(0 个)

类别

Help CenterFile Exchange 中查找有关 Lighting, Transparency, and Shading 的更多信息

Community Treasure Hunt

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

Start Hunting!

Translated by