How can I correct this code?
이전 댓글 표시
% Data Initialization
A = [ -0.1880
0.0100
-0.0079
-0.0273
0.1587
-0.4069
0.8943
-1.7404
3.0826
-5.1443
8.1501
-12.4334
18.3324
-26.3551
36.9611
-50.9219
68.8245
-91.8411
120.6744
-157.1298
201.9740
-258.0539
326.0616
-410.6317
512.0566
-638.1943
788.2017
-976.1147
1198.2528
-1481.4827
1815.1080
-2255.8411
2774.5508
-3510.8018
4379.5371
-5837.7578
7571.4097
-12936.2158
19555.4707
-10974.3643];
B = [ 11665.0576
-4.7293
13.9157
-26.5743
37.1578
-42.1387
44.1123
-45.5787
44.7030
-39.4265
30.5014
-20.7728
12.5115
-6.7371
3.2495
-1.4121
0.5494
-0.1905
0.0572
-0.0141
0.0023
0.0001
-0.0003
0.0002
-0.0001
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000];
C = [ -22.4193
-0.0217
0.0677
-0.1272
0.1340
-0.0548
-0.0433
0.1017
-0.1076
0.0619
-0.0047
-0.0195
0.0135
-0.0030
-0.0010
0.0009
-0.0002
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000];
D = [ -0.0871
-0.0006
0.0002
0.0004
-0.0012
0.0016
-0.0018
0.0017
-0.0015
0.0013
-0.0010
0.0008
-0.0006
0.0004
-0.0003
0.0002
-0.0001
0.0001
-0.0001
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000];
E = [ -0.1880
0.0100
-0.0079
-0.0273
0.1587
-0.4069
0.8943
-1.7404
3.0826
-5.1443
8.1501
-12.4334
18.3324
-26.3551
36.9611
-50.9219
68.8245
-91.8411
120.6744
-157.1298
201.9740
-258.0539
326.0616
-410.6317
512.0566
-638.1943
788.2017
-976.1147
1198.2528
-1481.4827
1815.1080
-2255.8411
2774.5508
-3510.8018
4379.5371
-5837.7578
7571.4097
-12936.2158
19555.4707
-10974.3643];
F = [ 11665.0576
-4.7293
13.9157
-26.5743
37.1578
-42.1387
44.1123
-45.5787
44.7030
-39.4265
30.5014
-20.7728
12.5115
-6.7371
3.2495
-1.4121
0.5494
-0.1905
0.0572
-0.0141
0.0023
0.0001
-0.0003
0.0002
-0.0001
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000];
% Parameters
a = 1; % Radius
L = 0.1; % Length
dd = 6; % Separation distance
alpha = 10;
etta = 0.01;
ettap = 0.1; % Additional parameters
zalm = alpha * complex(1, 0) / sqrt(2.d0);
xi = 1.0 / etta;
xi1 = 1.0 / ettap;
zk1 = sqrt((xi + sqrt(xi^2 - 4.d0 * xi * zalm^2)) / 2.0);
zk2 = sqrt((xi - sqrt(xi^2 - 4.d0 * xi * zalm^2)) / 2.0);
c = -a / L;
b = a / L;
m = a * 50; % Number of intervals
% Create Mesh Grid
[x, y] = meshgrid(c + dd : (b - c) / m : b, c : (b - c) / m : b);
% Apply Boundary Conditions
[I, J] = find(sqrt(x.^2 + y.^2) < (a - 0.01));
if ~isempty(I)
x(I, J) = NaN; % Use NaN to avoid affecting the visualization
y(I, J) = NaN;
end
% Compute Values
r = sqrt(x.^2 + y.^2);
t = atan2(y, x);
r2 = sqrt(r.^2 + dd.^2 - 2 .* r .* dd .* cos(t));
zet = (r.^2 - r2.^2 - dd.^2) ./ (2 .* r2 .* dd);
% Calculate Psi1
psi1 = zeros(size(x));
for i = 2:5
Ai = A(i-1);
Bi = B(i-1);
Ci = C(i-1);
Di = D(i-1);
Ei = E(i-1);
Fi = F(i-1);
psi1 = psi1 + (Ai .* r.^(-i + 1) + r.^(1 / 2) .* besselk(i - 1 / 2, r .* zk1) .* Bi + ...
r.^(1 / 2) .* besselk(i - 1 / 2, r .* zk2) .* Ci) .* ...
gegenbauerC(i, -1 / 2, cos(t)) + ...
(Di .* r2.^(-i + 1) + r2.^(1 / 2) .* besselk(i - 1 / 2, r2 .* zk1) .* Ei + ...
r2.^(1 / 2) .* besselk(i - 1 / 2, r2 .* zk2) .* Fi) .* ...
gegenbauerC(i, -1 / 2, zet);
end
% Plot Contours with Different Colors
figure;
contourf(y, x, psi1, 20); % Create filled contours with automatic color mapping
colorbar; % Add colorbar to check values
axis equal;
title('$\delta=1.5,\frac{\beta\mu}{a_1}=1.0,\alpha=1.0,\ell=0.1$', ...
'Interpreter', 'latex', 'FontSize', 12, 'FontName', 'Times New Roman', 'FontWeight', 'Normal');
% Plot Spheres
hold on;
% Sphere 1 (Vertical)
t3 = linspace(0, 2 * pi, 1000);
rr2 = 2.5;
x2 = rr2 * cos(t3); % Keep x and y aligned for vertical spheres
y2 = rr2 * sin(t3);
plot(y2, x2, '-k', 'LineWidth', 1.1);
fill(y2, x2, [0.5 0.5 0.5]); % Fill with gray color
% Sphere 2 (Vertical)
t2 = linspace(0, 2 * pi, 1000);
h = dd;
rr = 2.5;
x1 = rr * cos(t2) + h; % Keep x and y aligned for vertical spheres
y1 = rr * sin(t2);
plot(y1, x1, '-k', 'LineWidth', 1.1);
fill(y1, x1, [0.5 0.5 0.5]); % Fill with gray color
% Axis Formatting
axis off;
댓글 수: 7
Torsten
2025년 3월 19일
psi1 is complex-valued - you can't make a contour plot of complex-valued functions.
What you can plot is abs(psi1), real(psi1) and imag(psi1).
But you should first check whether psi1 is computed correctly in your code.
Fangjun Jiang
2025년 3월 20일
편집: Fangjun Jiang
2025년 3월 20일
See above. This line was missing. You can format the text to be code and run it right away to check.
dd = 6; % Separation distance
Cris LaPierre
2025년 3월 20일
Your updated code is missing the variable dd. Once I add that (copied from the code you originally posted), it doesn't error out, but there still is no contour plot. You code now results in psi1 being an array of NaNs.
Torsten
2025년 3월 21일
Try
[DH1,h1]=contour(y,x,real(psi1));
colorbar
and you will see that the range of your values for psi1 goes up until 1e154.
So a variation of the contour colors is restricted to a very small region in the graphics.
답변 (0개)
카테고리
도움말 센터 및 File Exchange에서 Line Plots에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!

