Clear Filters
Clear Filters

error in 3d plotting.

1 view (last 30 days)
tuhin
tuhin on 19 Mar 2024
Commented: tuhin on 1 Apr 2024
% Define the equations for r(x) and theta(x) soln
%fun = @(x, y, mu, lambda, ke, ko) [y(2); (-x * y(2) + y(1) - ke^2 * x^2 / (2 * mu + lambda) * y(1) + ko^2 * x^2 * y(3) / (2 * mu + lambda)) / x^2; ...
% y(4); (-x * y(4) + y(3) - ko^2 / mu * x^2 * y(1) - ke^2 / mu * x^2 * y(3)) / x^2];
% Define the equations as x*r(x) and x*theta(x) soln
fun = @(x, y, mu, lambda, ke, ko) [x * y(2);
(-x^2 * y(2) + x * y(1) - (ke^2 * x^4) / (2 * mu + lambda) * y(1) + (ko^2 * x^4 * y(3)) / (2 * mu + lambda)) / x^2;
x * y(4);
(-x^2 * y(4) + x * y(3) - (ko^2 * x^4) / mu * y(1) - (ke^2 * x^4) / mu * y(3)) / x^2];
% Define the initial conditions
initialConditions = [0.2; 0; 0; 0]; % r(0.75) = 0.1625, r'(0.75) = 0, Theta(0.75) = 0, Theta'(0.75) = 0
% Define the parameter ranges
mu_r_range = linspace(3, 20, 12);
lambda_r_range = linspace(50, 100, 10);
mu_theta_range = linspace(3, 20, 12);
lambda_theta_range = linspace(50, 100, 10);
ke_range = linspace(0.0001, 0.6, 10);
ko_range = linspace(0.0001, 0.6, 10);
% Initialize arrays to store peak extremum values
peak_r = zeros(length(ke_range), length(ko_range));
peak_theta = zeros(length(ke_range), length(ko_range));
% Solve the differential equations for each combination of parameters
for m = 1:length(ke_range)
for n = 1:length(ko_range)
ke = ke_range(m);
ko = ko_range(n);
% Initialize arrays to store extremum values for each solution
extremum_r = zeros(size(mu_r_range));
extremum_theta = zeros(size(mu_theta_range));
for i = 1:length(mu_r_range)
for j = 1:length(lambda_r_range)
for k = 1:length(mu_theta_range)
for l = 1:length(lambda_theta_range)
mu_r = mu_r_range(i);
lambda_r = lambda_r_range(j);
mu_theta = mu_theta_range(k);
lambda_theta = lambda_theta_range(l);
% Solve the differential equations
[x, y] = ode45(@(x, y) fun(x, y, mu_r, lambda_r, ke, ko), [0.75, 210], initialConditions);
% Calculate extremum values for r(x) and Theta(x)
extremum_r(i) = max(abs([max(y(:, 1)), min(y(:, 1))])); % Absolute maximum between max and min
extremum_theta(k) = max(abs([max(y(:, 3)), min(y(:, 3))])); % Absolute maximum between max and min
end
end
end
end
% Store the peak extremum values for each combination of ke and ko
peak_r(m, n) = max(extremum_r);
peak_theta(m, n) = max(extremum_theta);
end
end
% Plot the peak extremum values of r(x) for each combination of ke and ko
figure;
surf(ke_range, ko_range, peak_r);
xlabel('ke');
ylabel('ko');
zlabel('Peak Extremum Value of r(x)');
title('Peak Extremum Values of r(x)');
% Plot the peak extremum values of Theta(x) for each combination of ke and ko
figure;
surf(ke_range, ko_range, peak_theta);
xlabel('ke');
ylabel('ko');
zlabel('Peak Extremum Value of Theta(x)');
title('Peak Extremum Values of Theta(x)');
Want to 3d plot of x*r(x) and x*theta(x) with respect to ke and ko. But not getting any peak in the 3d plots as it was mentioned in the code. Any suggestions with possible errors?

Accepted Answer

Hassaan
Hassaan on 19 Mar 2024
% Define the ODE function according to your problem statement.
% It's important to ensure that this function correctly models the physical or mathematical system.
fun = @(x, y, mu, lambda, ke, ko) [x * y(2);
(-x^2 * y(2) + x * y(1) - (ke^2 * x^4) / (2 * mu + lambda) + (ko^2 * x^4 * y(3)) / (2 * mu + lambda)) / x^2;
x * y(4);
(-x^2 * y(4) + x * y(3) - (ko^2 * x^4) / mu * y(1) - (ke^2 * x^4) / mu * y(3)) / x^2];
% Define your initial conditions and ensure they are suitable for the problem.
initialConditions = [0.2; 0; 0; 0];
% Define parameter ranges. Ensure these ranges are logical for your problem.
ke_range = linspace(0.0001, 0.6, 10);
ko_range = linspace(0.0001, 0.6, 10);
% Preallocate matrices to store results for efficiency.
peak_r = zeros(length(ke_range), length(ko_range));
peak_theta = zeros(length(ke_range), length(ko_range));
% Loop through each combination of parameters.
for m = 1:length(ke_range)
for n = 1:length(ko_range)
ke = ke_range(m);
ko = ko_range(n);
% Define material properties or other parameters if they vary.
mu = 10; % Example value
lambda = 50; % Example value
% Solve the differential equations over a defined x domain.
% Adjust the x domain [0.75, 210] as necessary for your problem.
[x, y] = ode45(@(x, y) fun(x, y, mu, lambda, ke, ko), [0.75, 210], initialConditions);
% Extract and store the peak values for r(x) and theta(x).
% This method captures the peaks in the results for each ke, ko pair.
peak_r(m, n) = max(abs(y(:, 1)));
peak_theta(m, n) = max(abs(y(:, 3)));
end
end
% Visualization
% Create 3D surface plots to visualize how peak values change with ke and ko.
figure;
surf(ke_range, ko_range, peak_r');
title('Peak Extremum Values of r(x)');
xlabel('ke');
ylabel('ko');
zlabel('Peak r(x)');
colorbar;
figure;
surf(ke_range, ko_range, peak_theta');
title('Peak Extremum Values of Theta(x)');
xlabel('ke');
ylabel('ko');
zlabel('Peak Theta(x)');
colorbar;
-----------------------------------------------------------------------------------------------------------------------------------------------------
If you find the solution helpful and it resolves your issue, it would be greatly appreciated if you could accept the answer. Also, leaving an upvote and a comment are also wonderful ways to provide feedback.
It's important to note that the advice and code are based on limited information and meant for educational purposes. Users should verify and adapt the code to their specific needs, ensuring compatibility and adherence to ethical standards.
Professional Interests
  • Technical Services and Consulting
  • Embedded Systems | Firmware Developement | Simulations
  • Electrical and Electronics Engineering
Feel free to contact me.
  1 Comment
tuhin
tuhin on 1 Apr 2024
% Rearranged Data Matrix
dr_data = [2.1714, -0.0059213; 4.3429, 0.017543; 6.5143, 0.086; 8.6857, 0.15588; 10.8571, 0.20539; 13.0286, 0.19694; 15.2, 0.20442; 17.3714, 0.15659; 19.5429, 0.18918; 21.7143, 0.18262; 23.8857, 0.13818; 26.0571, 0.1539; 28.2286, 0.11195; 30.4, 0.12689; 32.5714, 0.068146; 34.7429, 0.028182; 36.9143, 0.013852; 39.0857, 0.039137; 41.2571, 0.00033664; 43.4286, -0.036782; 45.6, -0.043573; 47.7714, -0.060933; 49.9429, -0.030135; 52.1143, -0.043654; 54.2857, -0.039393; 56.4571, -0.030637; 58.6286, -0.03931; 60.8, -0.044883; 62.9714, -0.022349; 65.1429, -0.01046; 67.3143, 0.0014764; 69.4857, 0.012712; 71.6571, 0.0211; 73.8286, 0.02204];
dtheta_data = [2.1714, -0.011123; 4.3429, -0.31772; 6.5143, -0.3745; 8.6857, -0.40013; 10.8571, -0.50617; 13.0286, -0.49345; 15.2, -0.44292; 17.3714, -0.42858; 19.5429, -0.41354; 21.7143, -0.29636; 23.8857, -0.22671; 26.0571, -0.099143; 28.2286, 0.0087113; 30.4, -0.042737; 32.5714, 0.01474; 34.7429, 0.11353; 36.9143, 0.094084; 39.0857, 0.11759; 41.2571, 0.16252; 43.4286, 0.16718; 45.6, 0.22171; 47.7714, 0.25543; 49.9429, 0.25836; 52.1143, 0.23052; 54.2857, 0.12648; 56.4571, 0.15518; 58.6286, 0.20362; 60.8, 0.22967; 62.9714, 0.22242; 65.1429, 0.19043; 67.3143, 0.1749; 69.4857, 0.16304; 71.6571, 0.14256; 73.8286, 0.14299];
% Define parameters
lambda = 1.2;
R = 180; Or can use range anything above 75 to 400
rout = 75;
% Create figure for separate plots
for i = 1:length(kappa_range)
for j = 1:length(theta_k_range)
kappa = kappa_range(i);
theta_k = theta_k_range(j);
% Calculate functions for the chosen kappa and theta_k
omega_m = sqrt(kappa / (2 * (lambda + 1))) * sqrt((lambda + 2) * cos(theta_k) - sqrt((lambda + 2)^2 * cos(theta_k)^2 - 4 * (lambda + 1)));
omega_p = sqrt(kappa / (2 * (lambda + 1))) * sqrt((lambda + 2) * cos(theta_k) + sqrt((lambda + 2)^2 * cos(theta_k)^2 - 4 * (lambda + 1)));
A1 = (8 * R^2 * (kappa^2 - omega_m^2 * omega_p^2)) / ((rout^2 - 4 * R^2)^2 * (kappa^2 + omega_m^2 * omega_p^2)) * (-2 * (omega_m^2 + omega_p^2) / (omega_m^2 * omega_p^2) ...
+ rout * omega_m^2 * (kappa^2 - omega_p^4) / (kappa^2 * omega_p * (omega_m^2 - omega_p^2) * besselj(1, rout * omega_p)) ...
- rout * omega_p^2 * (kappa^2 - omega_m^4) / (kappa^2 * omega_m * (omega_m^2 - omega_p^2) * besselj(1, rout * omega_m)));
B1 = (8 * R^2 * (kappa^2 - omega_m^2 * omega_p^2) * sqrt((kappa^2 - omega_m^4) * (kappa^2 - omega_p^4))) / ((rout^2 - 4 * R^2)^2 * (kappa^2 + omega_m^2 * omega_p^2)) * ...
(2 / (omega_m^2 * omega_p^2 * kappa) + rout / (kappa * omega_m * (omega_m^2 - omega_p^2) * besselj(1, rout * omega_m)) ...
- rout / (kappa * omega_p * (omega_m^2 - omega_p^2) * besselj(1, rout * omega_p)));
% Define functions for fitting
dr = @(r, params) (2 * besselj(1, r * omega_p) ./ ((omega_m^2 - omega_p^2) * (kappa^2 + omega_m^2 * omega_p^2)) .* ...
((omega_m^2 * (kappa^2 - omega_p^4)) ./ (omega_p .* (params(1) + (16 * R^2 * (kappa^2 - omega_m^2 * omega_p^2)) ./ ...
(kappa^2 * omega_p^2 * (rout^2 - 4 * R^2)^2)) + (params(2) * omega_m^2 * omega_p .* sqrt((kappa^2 - omega_m^4) * (kappa^2 - omega_p^4))) ./ kappa))) ...
- (2 * besselj(1, r * omega_m) ./ ((omega_m^2 - omega_p^2) * (kappa^2 + omega_m^2 * omega_p^2)) .* ...
((omega_p^2 * (kappa^2 - omega_m^4)) ./ (omega_m .* (params(1) + (16 * R^2 * (kappa^2 - omega_m^2 * omega_p^2)) ./ ...
(kappa^2 * omega_m^2 * (rout^2 - 4 * R^2)^2)) + (params(2) * omega_m * omega_p^2 .* sqrt((kappa^2 - omega_m^4) * (kappa^2 - omega_p^4))) ./ kappa))) ...
- (16 * r * R^2 * (omega_m^2 + omega_p^2) .* (kappa^2 - omega_m^2 * omega_p^2)) ./ ...
(omega_p^2 * omega_m^2 * (rout^2 - 4 * R^2)^2 * (kappa^2 + omega_m^2 * omega_p^2));
dtheta = @(r, params) (2 * besselj(1, r * omega_p) ./ ((kappa^2 + omega_m^2 * omega_p^2) * (omega_m^2 - omega_p^2)) .* ...
(-sqrt((kappa^2 - omega_m^4) * (kappa^2 - omega_p^4)) ./ (omega_p .* (params(1) * kappa + (16 * R^2 * (kappa^2 - omega_m^2 * omega_p^2)) ./ ...
(kappa * omega_p^2 * (rout^2 - 4 * R^2)^2) - params(2) * omega_p * (kappa^2 - omega_m^4)))) ...
+ (2 * besselj(1, r * omega_m) ./ ((kappa^2 + omega_m^2 * omega_p^2) * (omega_m^2 - omega_p^2))) .* ...
(sqrt((kappa^2 - omega_m^4) * (kappa^2 - omega_p^4)) ./ (omega_m .* (params(1) * kappa + (16 * R^2 * (kappa^2 - omega_m^2 * omega_p^2)) ./ ...
(kappa * omega_m^2 * (rout^2 - 4 * R^2)^2) + params(2) * omega_m * (kappa^2 - omega_p^4)))) ...
+ (16 * r * R^2 * (kappa^2 - omega_m^2 * omega_p^2) * sqrt((kappa^2 - omega_m^4) * (kappa^2 - omega_p^4))) ./ ...
(kappa * omega_m^2 * omega_p^2 * (rout^2 - 4 * R^2)^2 * (kappa^2 + omega_m^2 * omega_p^2)));
% Fit dr data
p_dr = lsqcurvefit(@(params, r) dr(r, params), [omega_p, A1], dr_data(:, 1), dr_data(:, 2));
% Fit dtheta data
p_dtheta = lsqcurvefit(@(params, r) dtheta(r, params), [omega_p, A1], dtheta_data(:, 1), dtheta_data(:, 2));
% Plot dr and dtheta
figure;
subplot(2, 1, 1);
plot(r, dr(r, p_dr));
hold on;
scatter(dr_data(:, 1), dr_data(:, 2), 'r');
xlabel('r');
ylabel('dr');
title(['kappa = ', num2str(kappa), ', theta_k = ', num2str(theta_k)]);
legend('Fitted Curve', 'Simulation Data');
subplot(2, 1, 2);
plot(r, dtheta(r, p_dtheta));
hold on;
scatter(dtheta_data(:, 1), dtheta_data(:, 2), 'r');
xlabel('r');
ylabel('dtheta');
title(['kappa = ', num2str(kappa), ', theta_k = ', num2str(theta_k)]);
legend('Fitted Curve', 'Simulation Data');
end
end
Hi Hassaan,
I need your help with this technical problem. I want to fit the simulated data points of d_r(r) vs r and d_theta(r ) vs r using the below mentioned analytical solutions. There are two parameters. one is kappa and another one is theta_k. For the fitting the kappa can choose anything positive values the theta_k should be with in 0 to pi/2 any values. If required one can vary R also in between 75 to 1000 or so. That means lambda, kappa, theta_k, and R can be used as parameters to fit the datas. I am not getting any proper fitting of those two plots. I would be appreaciate any help or suggestion about ths fitting.

Sign in to comment.

More Answers (0)

Products


Release

R2023b

Community Treasure Hunt

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

Start Hunting!