Error in soln. Please help to solve the issue.

5 次查看（过去 30 天）
tuhin
tuhin 2024-9-7
评论： Walter Roberson 2024-9-7
% Define initial parameters
lambda_init = 1.36;
R_init = 500;
rout = 76.4;
% Define c
c = sqrt(4 - (rout/R_init)^2);
% Define k_e and k_o ranges
kappa_e_values = linspace(0.0, 0.05, 20);
kappa_o_values = linspace(0.0, 0.05, 20);
% Initialize arrays to store the results
max_rd_dr = zeros(length(kappa_e_values), length(kappa_o_values));
min_rd_dr = zeros(length(kappa_e_values), length(kappa_o_values));
amp_rd_dr = zeros(length(kappa_e_values), length(kappa_o_values));
mod_rd_dr = zeros(length(kappa_e_values), length(kappa_o_values));
max_rd_dtheta = zeros(length(kappa_e_values), length(kappa_o_values));
min_rd_dtheta = zeros(length(kappa_e_values), length(kappa_o_values));
amp_rd_dtheta = zeros(length(kappa_e_values), length(kappa_o_values));
mod_rd_dtheta = zeros(length(kappa_e_values), length(kappa_o_values));
% Loop over kappa_e and kappa_o values
for i = 1:length(kappa_e_values)
for j = 1:length(kappa_o_values)
kappa_e = kappa_e_values(i);
kappa_o = kappa_o_values(j);
% Calculate kappa and theta_k
kappa = sqrt(kappa_e^2 + kappa_o^2);
theta_k = atan(kappa_o / kappa_e);
% Initial guess for the solution
solinit = bvpinit(linspace(0.0001, rout, 100), [0, 0, 0, 0]);
% Solve the BVP
sol = bvp4c(@(r, y) odefun(r, y, lambda_init, kappa, theta_k, c, R_init), @bcfun, solinit);
% Extract solutions
r = linspace(0.0001, rout, 100);
y = deval(sol, r);
dr_sol = y(1,:);
dtheta_sol = y(3,:);
% Calculate r*dr and r*dtheta
r_dr = r .* dr_sol;
r_dtheta = r .* dtheta_sol;
% Calculate max, min, amplitude, and modulus of r*dr
max_rd_dr(i, j) = max(r_dr);
min_rd_dr(i, j) = min(r_dr);
amp_rd_dr(i, j) = max_rd_dr(i, j) - min_rd_dr(i, j);
mod_rd_dr(i, j) = max(abs(max_rd_dr(i, j)), abs(min_rd_dr(i, j)));
% Calculate max, min, amplitude, and modulus of r*dtheta
max_rd_dtheta(i, j) = max(r_dtheta);
min_rd_dtheta(i, j) = min(r_dtheta);
amp_rd_dtheta(i, j) = max_rd_dtheta(i, j) - min_rd_dtheta(i, j);
mod_rd_dtheta(i, j) = max(abs(max_rd_dtheta(i, j)), abs(min_rd_dtheta(i, j)));
end
end
% Define a custom colormap 'zet' that emphasizes peak values
zet_colormap = summer;
% Plot the results with kappa_e and kappa_o
subplot(1,2,1);
surf(kappa_e_values, kappa_o_values, amp_rd_dr');
xlabel('k_e');
ylabel('k_o');
zlabel('Amplitude of (r*d_r)');
title('Amplitude of r*d_r');
cb = colorbar;
colormap(zet_colormap);
caxis([min(amp_rd_dr(:)), max(amp_rd_dr(:))]);
subplot(1,2,2);
surf(kappa_e_values, kappa_o_values, amp_rd_dtheta');
xlabel('k_e');
ylabel('k_o');
zlabel('Amplitude of (r*d_\theta)');
title('Amplitude of r*d_\theta');
cb = colorbar;
colormap(zet_colormap);
caxis([min(amp_rd_dtheta(:)), max(amp_rd_dtheta(:))]);
% Define the function for the differential equations
function dydr = odefun(r, y, lambda_init, kappa, theta_k, c, R_init)
dydr = zeros(4,1);
dydr(1) = y(2);
dydr(2) = -((lambda_init+1)*y(2)+1/r*(kappa*r^2*cos(theta_k)-(lambda_init+1))*y(1)-kappa*r*y(3)*sin(theta_k)+16*lambda_init*r^2/(c^4*R_init^2))/(r*(lambda_init+1));
dydr(3) = y(4);
dydr(4) = -(y(4)+1/r*(kappa*r^2*cos(theta_k)-1)*y(3)+kappa*r*y(1)*sin(theta_k))/r;
end
% Boundary conditions
function res = bcfun(ya, yb)
res = [ya(1); ya(3); yb(1); yb(3)];
end
  2 个评论
Star Strider
Star Strider 2024-9-7
Looking at the ‘Extra_Parameters’ vector (that I added just before the bvp4c call):
Extra_Parameters = [lambda_init, kappa, theta_k, c, R_init]
‘theta_k’ is NaN, and displaying ‘dydr’ revealed its elements always to be either 0 or NaN.
Walter Roberson
Walter Roberson 2024-9-7
I suggest you replace
theta_k = atan(kappa_o / kappa_e);
with
theta_k = atan2(kappa_o, kappa_e);

请先登录，再进行评论。

请先登录，再回答此问题。

回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Programming 的更多信息

标签

产品


版本

R2023b

Community Treasure Hunt

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

Start Hunting!

Translated by