How can I correct this code?

조회 수: 2 (최근 30일)
Shreen El-Sapa
Shreen El-Sapa 2025년 3월 19일
댓글: Torsten 2025년 3월 21일
% 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
Error using contourf (line 55)
Input arguments must be real. Use the function REAL to get the real part of the inputs.
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
Shreen El-Sapa
Shreen El-Sapa 2025년 3월 20일
편집: Walter Roberson 2025년 3월 21일
clc
A =[ 0.4383
-0.4080
1.4130
-3.7497
8.7052
-18.3631
35.6095
-64.5347
110.4802
-180.8033
284.0856
-432.3704
638.2991
-920.6697
1296.9564
-1795.9568
2440.4854
-3274.5806
4326.2554
-5663.6147
7318.2793
-9397.8857
11932.9268
-15098.8994
18913.5957
-23675.0488
29361.4824
-36506.2852
44985.0234
-55821.1367
68631.0156
-85581.2891
105597.8906
-134030.5938
167689.9219
-224158.2500
291518.9688
-499382.0625
756808.3125
-425743.2188];
B=[ 784.3210
-41.3373
140.2780
-272.2328
368.2668
-379.0945
310.1994
-209.5842
118.3485
-54.6368
18.1556
-1.5761
-3.6361
3.8433
-2.6156
1.4482
-0.6952
0.2997
-0.1176
0.0427
-0.0144
0.0045
-0.0013
0.0004
-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];
C=[ 3.9761
-0.1226
0.3397
-0.5729
0.5937
-0.3249
0.0574
0.0449
-0.0379
0.0126
-0.0014
-0.0005
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
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000];
AA=[ -0.1536
0.0249
-0.0437
0.0587
-0.0688
0.0729
-0.0710
0.0645
-0.0553
0.0453
-0.0356
0.0271
-0.0200
0.0144
-0.0101
0.0070
-0.0048
0.0032
-0.0021
0.0014
-0.0009
0.0006
-0.0004
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];
BB=[ 6.3555
-0.2784
-0.1038
1.4075
-2.6351
2.8560
-2.3294
1.5033
-0.8018
0.3637
-0.1427
0.0493
-0.0152
0.0042
-0.0011
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
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000];
CC=[ -1.2158
0.1330
-0.2657
0.1976
-0.0445
-0.0175
0.0124
-0.0028
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
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000E+00
0.0000E+00
0.0000E+00];
a = 1 ; %RADIUS
L=.11;
etta=0.01; ettap=0.01; alpha=4;delta=1.5;
lambda= complex(0,alpha)/sqrt(2);
xi=sqrt(1./etta); xi1=sqrt(1./ettap);
alpha1=sqrt((xi+sqrt(xi^2-4.*xi.*lambda.^2))./2);
alpha2=sqrt((xi-sqrt(xi^2-4.*xi.*lambda.^2))./2);
dd=6; %h
c =-a/L;
b =a/L;
m =a*50; % NUMBER OF INTERVALS
[x,y]=meshgrid((c+dd:(b-c)/m:b),(c:(b-c)/m:b)');
[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
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);
warning on
psi1=0;
for i=2:25
Ai=A(i-1);Bi=B(i-1);Ci=C(i-1);AAi=AA(i-1);BBi=BB(i-1);CCi=CC(i-1);
%psi1=-psi1-(Ai.*r.^(-i-1)+r.^(-3./2).*besselk(i-1./2,r.*alpha1).*Bi+r.^(-3./2).*besselk(i-1./2,r.*alpha2).*Ci).*legendreP(i-1,cos(t))-(AAi.*r2.^(-i-1)+r2.^(-3./2).*besselk(i-1./2,r2.*alpha1).*BBi+r2.^(1./2).*besselk(i-1./2,r2.*alpha2).*CCi).*legendreP(i-1,zet);
psi1=psi1+(Ai.*r.^(-i+1)+r.^(1./2).*besselk(i-1./2,r.*alpha1).*Bi+r.^(1./2).*besselk(i-1./2,r.*alpha2).*Ci).*gegenbauerC(i,-1./2, cos(t))+(AAi.*r2.^(-i+1)+r2.^(1./2).*besselk(i-1./2,r2.*alpha1).*BBi+r2.^(1./2).*besselk(i-1./2,r2.*alpha2).*CCi).*gegenbauerC(i,-1./2,zet);
end
hold on
%[DH1,h1]=contour(x,y,real(psi1),'--k','LineWidth',1.1,'ShowText','on'); %,psi2,'--k',psi2,':k'
[DH1,h1]=contour(y,x,real(psi1),'-k');
%p1=contour(y,x,real(psi1),[0.001 0.001],'-c','LineWidth',1.3,'ShowText','on'); %,'ShowText','on'
%p2=contour(y,x,real(psi1),[0.005 0.005],'-r','LineWidth',1.3);
%p3=contour(y,x,real(psi1),[0.05 0.05],'-b','LineWidth',1.3);
%p4=contour(y,x,real(psi1),[.1 .1],'-k','LineWidth',1.3);
%p5=contour(y,x,real(psi1),[.2 .2],'--c','LineWidth',1.3);
%p6=contour(y,x,real(psi1),[.4 .4],'--r','LineWidth',1.3);
%p7=contour(y,x,real(psi1),[.6 .6],'--b','LineWidth',1.3);
%p8=contour(y,x,real(psi1),[.8 .8],'--k','LineWidth',1.3);
%p9=contour(y,x,real(psi1),[1 1],'-.c','LineWidth',1.3);
%p10=contour(y,x,real(psi1),[2 2],'-.r','LineWidth',1.3);
%p11=contour(y,x,real(psi1),[3 3],'-.b','LineWidth',1.3);
%p12=contour(y,x,real(psi1),[4 4],'-.k','LineWidth',1.3);
%axis equal; % Maintain aspect ratio
%axis tight; % Adjust axis limits to fit the data
%%%%%%%%%%%%%%%
hold on
t3 = linspace(0,2*pi,1000);
h2=0;
k2=0;
rr2=2;
x2 = rr2*cos(t3)+h2;
y2 = rr2*sin(t3)+k2;
set(plot(y2,x2,'-k'),'LineWidth',1.1);
fill(y2,x2,[0.7 0.7 0.7])
hold on
t2 = linspace(0,2*pi,1000);
h=dd;
k=0;
rr=1;
x1 = rr*cos(t2)+h;
y1 = rr*sin(t2)+k;
set(plot(y1,x1,'-k'),'LineWidth',1.1);
fill(y1,x1,[0.7 0.7 0.7])
axis off
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Torsten
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개)

카테고리

Help CenterFile Exchange에서 Line Plots에 대해 자세히 알아보기

Community Treasure Hunt

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

Start Hunting!

Translated by