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.
Shreen El-Sapa
2025 年 3 月 20 日
編集済み: Cris LaPierre
2025 年 3 月 20 日
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.
Shreen El-Sapa
2025 年 3 月 20 日
編集済み: Walter Roberson
2025 年 3 月 20 日
Shreen El-Sapa
2025 年 3 月 20 日
編集済み: Walter Roberson
2025 年 3 月 21 日
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!

