You must fit both data curves simultaneously, not in two separate calls to "lsqcurvefit".
Sie verfolgen jetzt diese Frage
- Aktualisierungen können Sie in Ihrem Feed verfolgter Inhalte sehen.
- Je nach Ihren Kommunikationseinstellungen können Sie auch E-Mails erhalten.
Simulation data Fittting problem
6 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
% 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
................................
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.
3 Kommentare
Torsten
am 1 Apr. 2024
⋮
-
Verknüpfen
Direkter Link zu diesem Kommentar
tuhin
am 1 Apr. 2024
⋮
-
Verknüpfen
Direkter Link zu diesem Kommentar
Thanks trosten! How to modify it? Can you please check the fitting with those 4 parmters. I checked and not working seems. If you have time please takae a look on the fitting.
tuhin
am 1 Apr. 2024
⋮
-
Verknüpfen
Direkter Link zu diesem Kommentar
Bearbeitet: Voss
am 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];
% Best Parameters (Assumed)
best_params = [0.2, 2, pi/10, 200]; % Just for illustration, replace with actual values
% Constants
lambda = best_params(1); % Best lambda value
kappa = best_params(2); % Best kappa value
theta_k = best_params(3); % Best theta_k value
R = best_params(4); % Best R value
rout = 75; % Define rout
% Calculate omega_m and omega_p
omega_m = @(lambda, kappa, theta_k) sqrt(kappa / (2 * (lambda + 1))) * sqrt((lambda + 2) * cos(theta_k) - sqrt((lambda + 2)^2 * cos(theta_k)^2 - 4 * (lambda + 1)));
omega_p = @(lambda, kappa, theta_k) sqrt(kappa / (2 * (lambda + 1))) * sqrt((lambda + 2) * cos(theta_k) + sqrt((lambda + 2)^2 * cos(theta_k)^2 - 4 * (lambda + 1)));
% Define A1 and B1
A1 = @(R, kappa, lambda, theta_k, rout) (8 * R^2 * (kappa^2 - omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2)) / ((rout^2 - 4 * R^2)^2 * (kappa^2 + omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2)) * ...
(-2 * (omega_m(lambda, kappa, theta_k)^2 + omega_p(lambda, kappa, theta_k)^2) / (omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2) ...
+ rout * omega_m(lambda, kappa, theta_k)^2 * (kappa^2 - omega_p(lambda, kappa, theta_k)^4) / (kappa^2 * omega_p(lambda, kappa, theta_k) * (omega_m(lambda, kappa, theta_k)^2 - omega_p(lambda, kappa, theta_k)^2) * besselj(1, rout * omega_p(lambda, kappa, theta_k))) ...
- rout * omega_p(lambda, kappa, theta_k)^2 * (kappa^2 - omega_m(lambda, kappa, theta_k)^4) / (kappa^2 * omega_m(lambda, kappa, theta_k) * (omega_m(lambda, kappa, theta_k)^2 - omega_p(lambda, kappa, theta_k)^2) * besselj(1, rout * omega_m(lambda, kappa, theta_k))));
B1 = @(R, kappa, lambda, theta_k, rout) (8 * R^2 * (kappa^2 - omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2) * sqrt((kappa^2 - omega_m(lambda, kappa, theta_k)^4) * (kappa^2 - omega_p(lambda, kappa, theta_k)^4))) / ((rout^2 - 4 * R^2)^2 * (kappa^2 + omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2)) * ...
(2 / (omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2 * kappa) + rout / (kappa * omega_m(lambda, kappa, theta_k) * (omega_m(lambda, kappa, theta_k)^2 - omega_p(lambda, kappa, theta_k)^2) * besselj(1, rout * omega_m(lambda, kappa, theta_k))) ...
- rout / (kappa * omega_p(lambda, kappa, theta_k) * (omega_m(lambda, kappa, theta_k)^2 - omega_p(lambda, kappa, theta_k)^2) * besselj(1, rout * omega_p(lambda, kappa, theta_k))));
% Define functions for fitting
dr = @(r, params) (2 * besselj(1, r * omega_p(lambda, kappa, theta_k)) ./ ((omega_m(lambda, kappa, theta_k)^2 - omega_p(lambda, kappa, theta_k)^2) * (kappa^2 + omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2)) .* ...
((omega_m(lambda, kappa, theta_k)^2 * (kappa^2 - omega_p(lambda, kappa, theta_k)^4)) ./ (omega_p(lambda, kappa, theta_k) .* (A1(R, kappa, lambda, theta_k, rout) + (16 * R^2 * (kappa^2 - omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2))) ./ ...
(kappa^2 * omega_p(lambda, kappa, theta_k)^2 * (rout^2 - 4 * R^2)^2) + (B1(R, kappa, lambda, theta_k, rout) * omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k) .* sqrt((kappa^2 - omega_m(lambda, kappa, theta_k)^4) * (kappa^2 - omega_p(lambda, kappa, theta_k)^4))) ./ kappa))) ...
- (2 * besselj(1, r * omega_m(lambda, kappa, theta_k)) ./ ((omega_m(lambda, kappa, theta_k)^2 - omega_p(lambda, kappa, theta_k)^2) * (kappa^2 + omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2)) .* ...
((omega_p(lambda, kappa, theta_k)^2 * (kappa^2 - omega_m(lambda, kappa, theta_k)^4)) ./ (omega_m(lambda, kappa, theta_k) .* (A1(R, kappa, lambda, theta_k, rout) + (16 * R^2 * (kappa^2 - omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2)) ./ ...
(kappa^2 * omega_m(lambda, kappa, theta_k)^2 * (rout^2 - 4 * R^2)^2) + (B1(R, kappa, lambda, theta_k, rout) * omega_m(lambda, kappa, theta_k) * omega_p(lambda, kappa, theta_k)^2 .* sqrt((kappa^2 - omega_m(lambda, kappa, theta_k)^4) * (kappa^2 - omega_p(lambda, kappa, theta_k)^4))) ./ kappa))) ...
- (16 * r * R^2 * (omega_m(lambda, kappa, theta_k)^2 + omega_p(lambda, kappa, theta_k)^2) .* (kappa^2 - omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2)) ./ ...
(omega_p(lambda, kappa, theta_k)^2 * omega_m(lambda, kappa, theta_k)^2 * (rout^2 - 4 * R^2)^2 * (kappa^2 + omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2)));
dtheta = @(r, params) (2 * besselj(1, r * omega_p(lambda, kappa, theta_k)) ./ ((kappa^2 + omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2) * (omega_m(lambda, kappa, theta_k)^2 - omega_p(lambda, kappa, theta_k)^2)) .* ...
(-sqrt((kappa^2 - omega_m(lambda, kappa, theta_k)^4) * (kappa^2 - omega_p(lambda, kappa, theta_k)^4)) ./ (omega_p(lambda, kappa, theta_k) .* (A1(R, kappa, lambda, theta_k, rout) * kappa + (16 * R^2 * (kappa^2 - omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2)) ./ ...
(kappa * omega_p(lambda, kappa, theta_k)^2 * (rout^2 - 4 * R^2)^2) - B1(R, kappa, lambda, theta_k, rout) * omega_p(lambda, kappa, theta_k) * (kappa^2 - omega_m(lambda, kappa, theta_k)^4)))) ...
+ (2 * besselj(1, r * omega_m(lambda, kappa, theta_k)) ./ ((kappa^2 + omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2) * (omega_m(lambda, kappa, theta_k)^2 - omega_p(lambda, kappa, theta_k)^2))) .* ...
(sqrt((kappa^2 - omega_m(lambda, kappa, theta_k)^4) * (kappa^2 - omega_p(lambda, kappa, theta_k)^4)) ./ (omega_m(lambda, kappa, theta_k) .* (A1(R, kappa, lambda, theta_k, rout) * kappa + (16 * R^2 * (kappa^2 - omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2)) ./ ...
(kappa * omega_m(lambda, kappa, theta_k)^2 * (rout^2 - 4 * R^2)^2) + B1(R, kappa, lambda, theta_k, rout) * omega_m(lambda, kappa, theta_k) * (kappa^2 - omega_p(lambda, kappa, theta_k)^4)))) ...
+ (16 * r * R^2 * (kappa^2 - omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2) * sqrt((kappa^2 - omega_m(lambda, kappa, theta_k)^4) * (kappa^2 - omega_p(lambda, kappa, theta_k)^4))) ./ ...
(kappa * omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2 * (rout^2 - 4 * R^2)^2 * (kappa^2 + omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2)));
% Fit the data using the defined functions with initial guesses from best_params
params_dr = lsqcurvefit(@(params, r) dr(r, params), best_params(1:2), dr_data(:, 1), dr_data(:, 2));
params_dtheta = lsqcurvefit(@(params, r) dtheta(r, params), best_params(1:2), dtheta_data(:, 1), dtheta_data(:, 2));
% Plot the results
figure;
subplot(1, 2, 1);
plot(dr_data(:, 1), dr_data(:, 2), 'ro', 'DisplayName', 'Data');
hold on;
plot(dr_data(:, 1), dr(dr_data(:, 1), params_dr), 'b-', 'DisplayName', 'Fit');
xlabel('r');
ylabel('dr');
legend('Location', 'best');
title('Fitting dr');
subplot(1, 2, 2);
plot(dtheta_data(:, 1), dtheta_data(:, 2), 'ro', 'DisplayName', 'Data');
hold on;
plot(dtheta_data(:, 1), dtheta(dtheta_data(:, 1), params_dtheta), 'b-', 'DisplayName', 'Fit');
xlabel('r');
ylabel('d\theta');
legend('Location', 'best');
title('Fitting d\theta');
Hi Torsten,Please check. still not working. please look into it
Antworten (1)
Torsten
am 1 Apr. 2024
Bearbeitet: Torsten
am 1 Apr. 2024
Call lsqcurvefit once as
dr = @(params,r) (2 * besselj(1, r(:,1). * 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(:,1) * 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(:,1) * 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 = @(params,r) (2 * besselj(1, r(:,2) * 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(:,2) * 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(:,2) * 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)));
params = lsqcurvefit(@(params,r)[dr(params,r),dtheta(params,r)],[omega_p, A1],[dr_data(:,1),dtheta_data(:,1)],[dr_data(:,2),dtheta_data(:,2)])
18 Kommentare
tuhin
am 1 Apr. 2024
⋮
-
Verknüpfen
Direkter Link zu diesem Kommentar
Can you please check the fitting once? Maybe there is some bug in the code. I want to get those parameter values from the fittings and want to use maxiterations for these.
Torsten
am 1 Apr. 2024
⋮
-
Verknüpfen
Direkter Link zu diesem Kommentar
Bearbeitet: Torsten
am 1 Apr. 2024
Matlab online does not work at the moment - thus no chance to run your code.
And please include the full code you are using at the moment with the changes from above incorporated.
tuhin
am 1 Apr. 2024
⋮
-
Verknüpfen
Direkter Link zu diesem Kommentar
% 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];
% Constants
lambda_init = best_params(1); % Best lambda value
kappa_init = best_params(2); % Best kappa value
theta_k_init = best_params(3); % Best theta_k value
R_init = best_params(4); % Best R value
rout = 75; % Define rout
% Calculate omega_m and omega_p
omega_m = @(lambda, kappa, theta_k) sqrt(kappa / (2 * (lambda + 1))) * sqrt((lambda + 2) * cos(theta_k) - sqrt((lambda + 2)^2 * cos(theta_k)^2 - 4 * (lambda + 1)));
omega_p = @(lambda, kappa, theta_k) sqrt(kappa / (2 * (lambda + 1))) * sqrt((lambda + 2) * cos(theta_k) + sqrt((lambda + 2)^2 * cos(theta_k)^2 - 4 * (lambda + 1)));
% Define A1 and B1
A1 = @(R, kappa, lambda, theta_k, rout) (8 * R^2 * (kappa^2 - omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2)) / ((rout^2 - 4 * R^2)^2 * (kappa^2 + omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2)) * ...
(-2 * (omega_m(lambda, kappa, theta_k)^2 + omega_p(lambda, kappa, theta_k)^2) / (omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2) ...
+ rout * omega_m(lambda, kappa, theta_k)^2 * (kappa^2 - omega_p(lambda, kappa, theta_k)^4) / (kappa^2 * omega_p(lambda, kappa, theta_k) * (omega_m(lambda, kappa, theta_k)^2 - omega_p(lambda, kappa, theta_k)^2) * besselj(1, rout * omega_p(lambda, kappa, theta_k))) ...
- rout * omega_p(lambda, kappa, theta_k)^2 * (kappa^2 - omega_m(lambda, kappa, theta_k)^4) / (kappa^2 * omega_m(lambda, kappa, theta_k) * (omega_m(lambda, kappa, theta_k)^2 - omega_p(lambda, kappa, theta_k)^2) * besselj(1, rout * omega_m(lambda, kappa, theta_k))));
B1 = @(R, kappa, lambda, theta_k, rout) (8 * R^2 * (kappa^2 - omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2) * sqrt((kappa^2 - omega_m(lambda, kappa, theta_k)^4) * (kappa^2 - omega_p(lambda, kappa, theta_k)^4))) / ((rout^2 - 4 * R^2)^2 * (kappa^2 + omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2) * ...
(kappa^2 + omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2));
% Specify the maximum number of iterations for fitting
max_iterations = 1000;
% Fit the data using lsqcurvefit with adjusted initial guesses and maximum iterations
options = optimoptions('lsqcurvefit', 'Algorithm', 'trust-region-reflective', 'MaxIterations', max_iterations);
% Define the anonymous functions for dr and dtheta
dr = @(r, params) (2 * besselj(1, r * omega_p(params(1), params(2), theta_k_init)) ./ ((params(2)^2 + omega_m(params(1), params(2), theta_k_init)^2 * omega_p(params(1), params(2), theta_k_init)^2) * (omega_m(params(1), params(2), theta_k_init)^2 - omega_p(params(1), params(2), theta_k_init)^2))) .* ...
((omega_m(params(1), params(2), theta_k_init)^2 * (params(2)^2 - omega_p(params(1), params(2), theta_k_init)^4)) ./ (omega_p(params(1), params(2), theta_k_init) .* (A1(R_init, params(2), params(1), theta_k_init, rout) * params(2) + (16 * R_init^2 * (params(2)^2 - omega_m(params(1), params(2), theta_k_init)^2 * omega_p(params(1), params(2), theta_k_init)^2)) ./ ...
(params(2) * omega_p(params(1), params(2), theta_k_init)^2 * (rout^2 - 4 * R_init^2)^2) - B1(R_init, params(2), params(1), theta_k_init, rout) * omega_p(params(1), params(2), theta_k_init) * (params(2)^2 - omega_m(params(1), params(2), theta_k_init)^4)))) ...
- (2 * besselj(1, r * omega_m(params(1), params(2), theta_k_init)) ./ ((params(2)^2 + omega_m(params(1), params(2), theta_k_init)^2 * omega_p(params(1), params(2), theta_k_init)^2) * (omega_m(params(1), params(2), theta_k_init)^2 - omega_p(params(1), params(2), theta_k_init)^2))) .* ...
((omega_p(params(1), params(2), theta_k_init)^2 * (params(2)^2 - omega_m(params(1), params(2), theta_k_init)^4)) ./ (omega_m(params(1), params(2), theta_k_init) .* (A1(R_init, params(2), params(1), theta_k_init, rout) * params(2) + (16 * R_init^2 * (params(2)^2 - omega_m(params(1), params(2), theta_k_init)^2 * omega_p(params(1), params(2), theta_k_init)^2)) ./ ...
(params(2) * omega_m(params(1), params(2), theta_k_init)^2 * (rout^2 - 4 * R_init^2)^2) + B1(R_init, params(2), params(1), theta_k_init, rout) * omega_m(params(1), params(2), theta_k_init) * (params(2)^2 - omega_p(params(1), params(2), theta_k_init)^4)))) ...
- (16 * r * R_init^2 * (omega_m(params(1), params(2), theta_k_init)^2 + omega_p(params(1), params(2), theta_k_init)^2) .* (params(2)^2 - omega_m(params(1), params(2), theta_k_init)^2 * omega_p(params(1), params(2), theta_k_init)^2)) ./ ...
(omega_p(params(1), params(2), theta_k_init)^2 * omega_m(params(1), params(2), theta_k_init)^2 * (rout^2 - 4 * R_init^2)^2 * (params(2)^2 + omega_m(params(1), params(2), theta_k_init)^2 * omega_p(params(1), params(2), theta_k_init)^2));
dtheta = @(r, params) (2 * besselj(1, r * omega_p(params(1), params(2), theta_k_init)) ./ ((params(2)^2 + omega_m(params(1), params(2), theta_k_init)^2 * omega_p(params(1), params(2), theta_k_init)^2) * (omega_m(params(1), params(2), theta_k_init)^2 - omega_p(params(1), params(2), theta_k_init)^2))) .* ...
(-sqrt((params(2)^2 - omega_m(params(1), params(2), theta_k_init)^4) * (params(2)^2 - omega_p(params(1), params(2), theta_k_init)^4)) ./ (omega_p(params(1), params(2), theta_k_init) .* (A1(R_init, params(2), params(1), theta_k_init, rout) * params(2) + (16 * R_init^2 * (params(2)^2 - omega_m(params(1), params(2), theta_k_init)^2 * omega_p(params(1), params(2), theta_k_init)^2)) ./ ...
(params(2) * omega_p(params(1), params(2), theta_k_init)^2 * (rout^2 - 4 * R_init^2)^2) - B1(R_init, params(2), params(1), theta_k_init, rout) * omega_p(params(1), params(2), theta_k_init) * (params(2)^2 - omega_m(params(1), params(2), theta_k_init)^4)))) ...
+ (2 * besselj(1, r * omega_m(params(1), params(2), theta_k_init)) ./ ((params(2)^2 + omega_m(params(1), params(2), theta_k_init)^2 * omega_p(params(1), params(2), theta_k_init)^2) * (omega_m(params(1), params(2), theta_k_init)^2 - omega_p(params(1), params(2), theta_k_init)^2))) .* ...
(sqrt((params(2)^2 - omega_m(params(1), params(2), theta_k_init)^4) * (params(2)^2 - omega_p(params(1), params(2), theta_k_init)^4)) ./ (omega_m(params(1), params(2), theta_k_init) .* (A1(R_init, params(2), params(1), theta_k_init, rout) * params(2) + (16 * R_init^2 * (params(2)^2 - omega_m(params(1), params(2), theta_k_init)^2 * omega_p(params(1), params(2), theta_k_init)^2)) ./ ...
(params(2) * omega_m(params(1), params(2), theta_k_init)^2 * (rout^2 - 4 * R_init^2)^2) + B1(R_init, params(2), params(1), theta_k_init, rout) * omega_m(params(1), params(2), theta_k_init) * (params(2)^2 - omega_p(params(1), params(2), theta_k_init)^4)))) ...
+ (16 * r * R_init^2 * (params(2)^2 - omega_m(params(1), params(2), theta_k_init)^2 * omega_p(params(1), params(2), theta_k_init)^2) * sqrt((params(2)^2 - omega_m(params(1), params(2), theta_k_init)^4) * (params(2)^2 - omega_p(params(1), params(2), theta_k_init)^4))) ./ ...
(params(2) * omega_m(params(1), params(2), theta_k_init)^2 * omega_p(params(1), params(2), theta_k_init)^2 * (rout^2 - 4 * R_init^2)^2 * (params(2)^2 + omega_m(params(1), params(2), theta_k_init)^2 * omega_p(params(1), params(2), theta_k_init)^2));
% Fit the data using the defined functions with initial guesses from best_params and specified options
params = lsqcurvefit(@(params, r) [dr(r, params); dtheta(r, params)], [lambda_init, kappa_init], [dr_data(:, 1); dtheta_data(:, 1)], [dr_data(:, 2); dtheta_data(:, 2)], [], [], options);
% Extract the optimized parameters
lambda_fit = params(1);
kappa_fit = params(2);
theta_k_fit = theta_k_init; % Use the initial guess for theta_k
R_fit = R_init;
% Plotting data and fitting
r_values = linspace(min([dr_data(:, 1); dtheta_data(:, 1)]), max([dr_data(:, 1); dtheta_data(:, 1)]), 1000);
dr_fit = dr(r_values, params);
dtheta_fit = dtheta(r_values, params);
figure;
subplot(2, 1, 1);
plot(dr_data(:, 1), dr_data(:, 2), 'bo', 'DisplayName', 'Data');
hold on;
plot(r_values, dr_fit, 'r-', 'DisplayName', 'Fitted');
xlabel('r');
ylabel('dr');
title('Fitting dr');
legend;
subplot(2, 1, 2);
plot(dtheta_data(:, 1), dtheta_data(:, 2), 'bo', 'DisplayName', 'Data');
hold on;
plot(r_values, dtheta_fit, 'r-', 'DisplayName', 'Fitted');
xlabel('r');
ylabel('dtheta');
title('Fitting dtheta');
legend;
getting error while run it. but as you suggested i modified it. I am expecting to the analytic solutions to fit the data and from fitting Iwill get those paramaeter vlues. Please check it once for me. It would be helpful
Torsten
am 1 Apr. 2024
⋮
-
Verknüpfen
Direkter Link zu diesem Kommentar
Bearbeitet: Torsten
am 1 Apr. 2024
Look into my code about how r is referenced in the definition of "dr" and "dtheta". It's wrong in your code.
And "best_params" is undefined at the beginning.
Check whether your function works properly by calling it for the initial parameter values as
init = [dr([dr_data(:, 1); dtheta_data(:, 1)],[lambda_init, kappa_init]);dtheta([dr_data(:, 1); dtheta_data(:, 1)],[lambda_init, kappa_init])]
before calling "lsqcurvefit".
tuhin
am 1 Apr. 2024
⋮
-
Verknüpfen
Direkter Link zu diesem Kommentar
I didnot get you. cna you please write the fullone.
However, this one seems working butnot fitting properly due to the proper parameter values variataions.
here is the one.you can look and check this one:
% 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];
% Manually specified initial guess values for parameters
lambda_init = 0.5; % Example initial guess for lambda
kappa_init = 0.5; % Example initial guess for kappa
theta_k_init = 0.1; % Example initial guess for theta_k
R_init = 50; % Example initial guess for R
% Constants
rout = 75; % Define rout
% Calculate omega_m and omega_p
omega_m = @(lambda, kappa, theta_k) sqrt(kappa / (2 * (lambda + 1))) * sqrt((lambda + 2) * cos(theta_k) - sqrt((lambda + 2)^2 * cos(theta_k)^2 - 4 * (lambda + 1)));
omega_p = @(lambda, kappa, theta_k) sqrt(kappa / (2 * (lambda + 1))) * sqrt((lambda + 2) * cos(theta_k) + sqrt((lambda + 2)^2 * cos(theta_k)^2 - 4 * (lambda + 1)));
% Define A1 and B1
A1 = @(R, kappa, lambda, theta_k, rout) (8 * R^2 * (kappa^2 - omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2)) / ((rout^2 - 4 * R^2)^2 * (kappa^2 + omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2)) * ...
(-2 * (omega_m(lambda, kappa, theta_k)^2 + omega_p(lambda, kappa, theta_k)^2) / (omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2) ...
+ rout * omega_m(lambda, kappa, theta_k)^2 * (kappa^2 - omega_p(lambda, kappa, theta_k)^4) / (kappa^2 * omega_p(lambda, kappa, theta_k) * (omega_m(lambda, kappa, theta_k)^2 - omega_p(lambda, kappa, theta_k)^2) * besselj(1, rout * omega_p(lambda, kappa, theta_k))) ...
- rout * omega_p(lambda, kappa, theta_k)^2 * (kappa^2 - omega_m(lambda, kappa, theta_k)^4) / (kappa^2 * omega_m(lambda, kappa, theta_k) * (omega_m(lambda, kappa, theta_k)^2 - omega_p(lambda, kappa, theta_k)^2) * besselj(1, rout * omega_m(lambda, kappa, theta_k))));
B1 = @(R, kappa, lambda, theta_k, rout) (8 * R^2 * (kappa^2 - omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2) * sqrt((kappa^2 - omega_m(lambda, kappa, theta_k)^4) * (kappa^2 - omega_p(lambda, kappa, theta_k)^4))) / ((rout^2 - 4 * R^2)^2 * (kappa^2 + omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2) * ...
(kappa^2 + omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2));
dr = @(r, params) (2 * besselj(1, r * omega_p(lambda, kappa, theta_k)) ./ ((omega_m(lambda, kappa, theta_k)^2 - omega_p(lambda, kappa, theta_k)^2) * (kappa^2 + omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2)) .* ...
((omega_m(lambda, kappa, theta_k)^2 * (kappa^2 - omega_p(lambda, kappa, theta_k)^4)) ./ (omega_p(lambda, kappa, theta_k) .* (A1(R, kappa, lambda, theta_k, rout) + (16 * R^2 * (kappa^2 - omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2))) ./ ...
(kappa^2 * omega_p(lambda, kappa, theta_k)^2 * (rout^2 - 4 * R^2)^2) + (B1(R, kappa, lambda, theta_k, rout) * omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k) .* sqrt((kappa^2 - omega_m(lambda, kappa, theta_k)^4) * (kappa^2 - omega_p(lambda, kappa, theta_k)^4))) ./ kappa))) ...
- (2 * besselj(1, r * omega_m(lambda, kappa, theta_k)) ./ ((omega_m(lambda, kappa, theta_k)^2 - omega_p(lambda, kappa, theta_k)^2) * (kappa^2 + omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2)) .* ...
((omega_p(lambda, kappa, theta_k)^2 * (kappa^2 - omega_m(lambda, kappa, theta_k)^4)) ./ (omega_m(lambda, kappa, theta_k) .* (A1(R, kappa, lambda, theta_k, rout) + (16 * R^2 * (kappa^2 - omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2)) ./ ...
(kappa^2 * omega_m(lambda, kappa, theta_k)^2 * (rout^2 - 4 * R^2)^2) + (B1(R, kappa, lambda, theta_k, rout) * omega_m(lambda, kappa, theta_k) * omega_p(lambda, kappa, theta_k)^2 .* sqrt((kappa^2 - omega_m(lambda, kappa, theta_k)^4) * (kappa^2 - omega_p(lambda, kappa, theta_k)^4))) ./ kappa))) ...
- (16 * r * R^2 * (omega_m(lambda, kappa, theta_k)^2 + omega_p(lambda, kappa, theta_k)^2) .* (kappa^2 - omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2)) ./ ...
(omega_p(lambda, kappa, theta_k)^2 * omega_m(lambda, kappa, theta_k)^2 * (rout^2 - 4 * R^2)^2 * (kappa^2 + omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2)));
dtheta = @(r, params) (2 * besselj(1, r * omega_p(lambda, kappa, theta_k)) ./ ((kappa^2 + omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2) * (omega_m(lambda, kappa, theta_k)^2 - omega_p(lambda, kappa, theta_k)^2)) .* ...
(-sqrt((kappa^2 - omega_m(lambda, kappa, theta_k)^4) * (kappa^2 - omega_p(lambda, kappa, theta_k)^4)) ./ (omega_p(lambda, kappa, theta_k) .* (A1(R, kappa, lambda, theta_k, rout) * kappa + (16 * R^2 * (kappa^2 - omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2)) ./ ...
(kappa * omega_p(lambda, kappa, theta_k)^2 * (rout^2 - 4 * R^2)^2) - B1(R, kappa, lambda, theta_k, rout) * omega_p(lambda, kappa, theta_k) * (kappa^2 - omega_m(lambda, kappa, theta_k)^4)))) ...
+ (2 * besselj(1, r * omega_m(lambda, kappa, theta_k)) ./ ((kappa^2 + omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2) * (omega_m(lambda, kappa, theta_k)^2 - omega_p(lambda, kappa, theta_k)^2))) .* ...
(sqrt((kappa^2 - omega_m(lambda, kappa, theta_k)^4) * (kappa^2 - omega_p(lambda, kappa, theta_k)^4)) ./ (omega_m(lambda, kappa, theta_k) .* (A1(R, kappa, lambda, theta_k, rout) * kappa + (16 * R^2 * (kappa^2 - omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2)) ./ ...
(kappa * omega_m(lambda, kappa, theta_k)^2 * (rout^2 - 4 * R^2)^2) + B1(R, kappa, lambda, theta_k, rout) * omega_m(lambda, kappa, theta_k) * (kappa^2 - omega_p(lambda, kappa, theta_k)^4)))) ...
+ (16 * r * R^2 * (kappa^2 - omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2) * sqrt((kappa^2 - omega_m(lambda, kappa, theta_k)^4) * (kappa^2 - omega_p(lambda, kappa, theta_k)^4))) ./ ...
(kappa * omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2 * (rout^2 - 4 * R^2)^2 * (kappa^2 + omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2)));
% Specify the maximum number of iterations for fitting
max_iterations = 1000000;
% Fit the data for dr using lsqcurvefit
options_dr = optimoptions('lsqcurvefit', 'Algorithm', 'trust-region-reflective', 'MaxIterations', max_iterations);
params_dr = lsqcurvefit(@(params, r) [dr(r, params)], [lambda_init, kappa_init], dr_data(:, 1), dr_data(:, 2), [], [], options_dr);
% Fit the data for dtheta using lsqcurvefit
options_dtheta = optimoptions('lsqcurvefit', 'Algorithm', 'trust-region-reflective', 'MaxIterations', max_iterations);
params_dtheta = lsqcurvefit(@(params, r) [dtheta(r, params)], [lambda_init, kappa_init], dtheta_data(:, 1), dtheta_data(:, 2), [], [], options_dtheta);
% Extract the optimized parameters for dr
lambda_fit_dr = params_dr(1);
kappa_fit_dr = params_dr(2);
% Extract the optimized parameters for dtheta
lambda_fit_dtheta = params_dtheta(1);
kappa_fit_dtheta = params_dtheta(2);
% Use the initial guess for theta_k and R
theta_k_fit = theta_k_init;
R_fit = R_init;
% Plotting data and fitting for dr
r_values_dr = linspace(min(dr_data(:, 1)), max(dr_data(:, 1)), 1000);
dr_fit = dr(r_values_dr, params_dr);
% Plotting data and fitting for dtheta
r_values_dtheta = linspace(min(dtheta_data(:, 1)), max(dtheta_data(:, 1)), 1000);
dtheta_fit = dtheta(r_values_dtheta, params_dtheta);
% Plotting
figure;
subplot(2, 1, 1);
plot(dr_data(:, 1), dr_data(:, 2), 'bo', 'DisplayName', 'Data');
hold on;
plot(r_values_dr, dr_fit, 'r-', 'DisplayName', 'Fitted');
xlabel('r');
ylabel('dr');
title('Fitting dr');
legend;
subplot(2, 1, 2);
plot(dtheta_data(:, 1), dtheta_data(:, 2), 'bo', 'DisplayName', 'Data');
hold on;
plot(r_values_dtheta, dtheta_fit, 'r-', 'DisplayName', 'Fitted');
xlabel('r');
ylabel('dtheta');
title('Fitting dtheta');
legend;
tuhin
am 1 Apr. 2024
⋮
-
Verknüpfen
Direkter Link zu diesem Kommentar
Bearbeitet: tuhin
am 1 Apr. 2024
For your reference the anatic formula is given below for dr and dtheta. you can take a look into it. The analytic soln mentioned in the pdf is written in the aabove codenad trying to use this to fit the data. Please check this once.
Torsten
am 1 Apr. 2024
⋮
-
Verknüpfen
Direkter Link zu diesem Kommentar
Bearbeitet: Torsten
am 1 Apr. 2024
Don't you understand that you have to get the same parameters in the calls to dr and dtheta and that you consequently have to fit both functions simultaneously (in one call to lsqcurvefit) as shown in my code fragment ?
tuhin
am 1 Apr. 2024
⋮
-
Verknüpfen
Direkter Link zu diesem Kommentar
yes, but that not working again. here is the updated code with one call:
% 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];
% Constants
best_params = [1.2, 1,pi/4,150];
lambda_init = best_params(1); % Best lambda value
kappa_init = best_params(2); % Best kappa value
theta_k_init = best_params(3); % Best theta_k value
R_init = best_params(4); % Best R value
rout = 75; % Define rout
% Calculate omega_m and omega_p
omega_m = @(lambda, kappa, theta_k) sqrt(kappa / (2 * (lambda + 1))) * sqrt((lambda + 2) * cos(theta_k) - sqrt((lambda + 2)^2 * cos(theta_k)^2 - 4 * (lambda + 1)));
omega_p = @(lambda, kappa, theta_k) sqrt(kappa / (2 * (lambda + 1))) * sqrt((lambda + 2) * cos(theta_k) + sqrt((lambda + 2)^2 * cos(theta_k)^2 - 4 * (lambda + 1)));
% Define A1 and B1
A1 = @(R, kappa, lambda, theta_k, rout) (8 * R^2 * (kappa^2 - omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2)) / ((rout^2 - 4 * R^2)^2 * (kappa^2 + omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2)) * ...
(-2 * (omega_m(lambda, kappa, theta_k)^2 + omega_p(lambda, kappa, theta_k)^2) / (omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2) ...
+ rout * omega_m(lambda, kappa, theta_k)^2 * (kappa^2 - omega_p(lambda, kappa, theta_k)^4) / (kappa^2 * omega_p(lambda, kappa, theta_k) * (omega_m(lambda, kappa, theta_k)^2 - omega_p(lambda, kappa, theta_k)^2) * besselj(1, rout * omega_p(lambda, kappa, theta_k))) ...
- rout * omega_p(lambda, kappa, theta_k)^2 * (kappa^2 - omega_m(lambda, kappa, theta_k)^4) / (kappa^2 * omega_m(lambda, kappa, theta_k) * (omega_m(lambda, kappa, theta_k)^2 - omega_p(lambda, kappa, theta_k)^2) * besselj(1, rout * omega_m(lambda, kappa, theta_k))));
B1 = @(R, kappa, lambda, theta_k, rout) (8 * R^2 * (kappa^2 - omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2) * sqrt((kappa^2 - omega_m(lambda, kappa, theta_k)^4) * (kappa^2 - omega_p(lambda, kappa, theta_k)^4))) / ((rout^2 - 4 * R^2)^2 * (kappa^2 + omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2) * ...
(kappa^2 + omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2));
% Specify the maximum number of iterations for fitting
max_iterations = 1000;
% Fit the data using lsqcurvefit with adjusted initial guesses and maximum iterations
options = optimoptions('lsqcurvefit', 'Algorithm', 'trust-region-reflective', 'MaxIterations', max_iterations);
% Define the anonymous functions for dr and dtheta
dr = @(r, params) (2 * besselj(1, r * omega_p(params(1), params(2), theta_k_init)) ./ ((params(2)^2 + omega_m(params(1), params(2), theta_k_init)^2 * omega_p(params(1), params(2), theta_k_init)^2) * (omega_m(params(1), params(2), theta_k_init)^2 - omega_p(params(1), params(2), theta_k_init)^2))) .* ...
((omega_m(params(1), params(2), theta_k_init)^2 * (params(2)^2 - omega_p(params(1), params(2), theta_k_init)^4)) ./ (omega_p(params(1), params(2), theta_k_init) .* (A1(R_init, params(2), params(1), theta_k_init, rout) * params(2) + (16 * R_init^2 * (params(2)^2 - omega_m(params(1), params(2), theta_k_init)^2 * omega_p(params(1), params(2), theta_k_init)^2)) ./ ...
(params(2) * omega_p(params(1), params(2), theta_k_init)^2 * (rout^2 - 4 * R_init^2)^2) - B1(R_init, params(2), params(1), theta_k_init, rout) * omega_p(params(1), params(2), theta_k_init) * (params(2)^2 - omega_m(params(1), params(2), theta_k_init)^4)))) ...
- (2 * besselj(1, r * omega_m(params(1), params(2), theta_k_init)) ./ ((params(2)^2 + omega_m(params(1), params(2), theta_k_init)^2 * omega_p(params(1), params(2), theta_k_init)^2) * (omega_m(params(1), params(2), theta_k_init)^2 - omega_p(params(1), params(2), theta_k_init)^2))) .* ...
((omega_p(params(1), params(2), theta_k_init)^2 * (params(2)^2 - omega_m(params(1), params(2), theta_k_init)^4)) ./ (omega_m(params(1), params(2), theta_k_init) .* (A1(R_init, params(2), params(1), theta_k_init, rout) * params(2) + (16 * R_init^2 * (params(2)^2 - omega_m(params(1), params(2), theta_k_init)^2 * omega_p(params(1), params(2), theta_k_init)^2)) ./ ...
(params(2) * omega_m(params(1), params(2), theta_k_init)^2 * (rout^2 - 4 * R_init^2)^2) + B1(R_init, params(2), params(1), theta_k_init, rout) * omega_m(params(1), params(2), theta_k_init) * (params(2)^2 - omega_p(params(1), params(2), theta_k_init)^4)))) ...
- (16 * r * R_init^2 * (omega_m(params(1), params(2), theta_k_init)^2 + omega_p(params(1), params(2), theta_k_init)^2) .* (params(2)^2 - omega_m(params(1), params(2), theta_k_init)^2 * omega_p(params(1), params(2), theta_k_init)^2)) ./ ...
(omega_p(params(1), params(2), theta_k_init)^2 * omega_m(params(1), params(2), theta_k_init)^2 * (rout^2 - 4 * R_init^2)^2 * (params(2)^2 + omega_m(params(1), params(2), theta_k_init)^2 * omega_p(params(1), params(2), theta_k_init)^2));
dtheta = @(r, params) (2 * besselj(1, r * omega_p(params(1), params(2), theta_k_init)) ./ ((params(2)^2 + omega_m(params(1), params(2), theta_k_init)^2 * omega_p(params(1), params(2), theta_k_init)^2) * (omega_m(params(1), params(2), theta_k_init)^2 - omega_p(params(1), params(2), theta_k_init)^2))) .* ...
(-sqrt((params(2)^2 - omega_m(params(1), params(2), theta_k_init)^4) * (params(2)^2 - omega_p(params(1), params(2), theta_k_init)^4)) ./ (omega_p(params(1), params(2), theta_k_init) .* (A1(R_init, params(2), params(1), theta_k_init, rout) * params(2) + (16 * R_init^2 * (params(2)^2 - omega_m(params(1), params(2), theta_k_init)^2 * omega_p(params(1), params(2), theta_k_init)^2)) ./ ...
(params(2) * omega_p(params(1), params(2), theta_k_init)^2 * (rout^2 - 4 * R_init^2)^2) - B1(R_init, params(2), params(1), theta_k_init, rout) * omega_p(params(1), params(2), theta_k_init) * (params(2)^2 - omega_m(params(1), params(2), theta_k_init)^4)))) ...
+ (2 * besselj(1, r * omega_m(params(1), params(2), theta_k_init)) ./ ((params(2)^2 + omega_m(params(1), params(2), theta_k_init)^2 * omega_p(params(1), params(2), theta_k_init)^2) * (omega_m(params(1), params(2), theta_k_init)^2 - omega_p(params(1), params(2), theta_k_init)^2))) .* ...
(sqrt((params(2)^2 - omega_m(params(1), params(2), theta_k_init)^4) * (params(2)^2 - omega_p(params(1), params(2), theta_k_init)^4)) ./ (omega_m(params(1), params(2), theta_k_init) .* (A1(R_init, params(2), params(1), theta_k_init, rout) * params(2) + (16 * R_init^2 * (params(2)^2 - omega_m(params(1), params(2), theta_k_init)^2 * omega_p(params(1), params(2), theta_k_init)^2)) ./ ...
(params(2) * omega_m(params(1), params(2), theta_k_init)^2 * (rout^2 - 4 * R_init^2)^2) + B1(R_init, params(2), params(1), theta_k_init, rout) * omega_m(params(1), params(2), theta_k_init) * (params(2)^2 - omega_p(params(1), params(2), theta_k_init)^4)))) ...
+ (16 * r * R_init^2 * (params(2)^2 - omega_m(params(1), params(2), theta_k_init)^2 * omega_p(params(1), params(2), theta_k_init)^2) * sqrt((params(2)^2 - omega_m(params(1), params(2), theta_k_init)^4) * (params(2)^2 - omega_p(params(1), params(2), theta_k_init)^4))) ./ ...
(params(2) * omega_m(params(1), params(2), theta_k_init)^2 * omega_p(params(1), params(2), theta_k_init)^2 * (rout^2 - 4 * R_init^2)^2 * (params(2)^2 + omega_m(params(1), params(2), theta_k_init)^2 * omega_p(params(1), params(2), theta_k_init)^2));
% Fit the data using the defined functions with initial guesses from best_params and specified options
init = [dr([dr_data(:, 1); dtheta_data(:, 1)], [lambda_init, kappa_init]); dtheta([dr_data(:, 1); dtheta_data(:, 1)], [lambda_init, kappa_init])];
params = lsqcurvefit(@(params, r) [dr(r, params); dtheta(r, params)], [lambda_init, kappa_init], [dr_data(:, 1); dtheta_data(:, 1)], [dr_data(:, 2); dtheta_data(:, 2)], [], [], options);
% Extract the optimized parameters
lambda_fit = params(1);
kappa_fit = params(2);
theta_k_fit = theta_k_init; % Use the initial guess for theta_k
R_fit = R_init;
% Plotting data and fitting
r_values = linspace(min([dr_data(:, 1); dtheta_data(:, 1)]), max([dr_data(:, 1); dtheta_data(:, 1)]), 1000);
dr_fit = dr(r_values, params);
dtheta_fit = dtheta(r_values, params);
figure;
subplot(2, 1, 1);
plot(dr_data(:, 1), dr_data(:, 2), 'bo', 'DisplayName', 'Data');
hold on;
plot(r_values, dr_fit, 'r-', 'DisplayName', 'Fitted');
xlabel('r');
ylabel('dr');
title('Fitting dr');
legend;
subplot(2, 1, 2);
plot(dtheta_data(:, 1), dtheta_data(:, 2), 'bo', 'DisplayName', 'Data');
hold on;
plot(r_values, dtheta_fit, 'r-', 'DisplayName', 'Fitted');
xlabel('r');
ylabel('dtheta');
title('Fitting dtheta');
legend;
Torsten
am 2 Apr. 2024
⋮
-
Verknüpfen
Direkter Link zu diesem Kommentar
Bearbeitet: Torsten
am 2 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];
% Manually specified initial guess values for parameters
lambda_init = 0.5; % Example initial guess for lambda
kappa_init = 0.5; % Example initial guess for kappa
theta_k = 0.1; % Example initial guess for theta_k
R = 50; % Example initial guess for R
% Constants
rout = 75; % Define rout
% Calculate omega_m and omega_p
omega_m = @(lambda, kappa, theta_k) sqrt(kappa / (2 * (lambda + 1))) * sqrt((lambda + 2) * cos(theta_k) - sqrt((lambda + 2)^2 * cos(theta_k)^2 - 4 * (lambda + 1)));
omega_p = @(lambda, kappa, theta_k) sqrt(kappa / (2 * (lambda + 1))) * sqrt((lambda + 2) * cos(theta_k) + sqrt((lambda + 2)^2 * cos(theta_k)^2 - 4 * (lambda + 1)));
% Define A1 and B1
A1 = @(R, kappa, lambda, theta_k, rout) (8 * R^2 * (kappa^2 - omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2)) / ((rout^2 - 4 * R^2)^2 * (kappa^2 + omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2)) * ...
(-2 * (omega_m(lambda, kappa, theta_k)^2 + omega_p(lambda, kappa, theta_k)^2) / (omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2) ...
+ rout * omega_m(lambda, kappa, theta_k)^2 * (kappa^2 - omega_p(lambda, kappa, theta_k)^4) / (kappa^2 * omega_p(lambda, kappa, theta_k) * (omega_m(lambda, kappa, theta_k)^2 - omega_p(lambda, kappa, theta_k)^2) * besselj(1, rout * omega_p(lambda, kappa, theta_k))) ...
- rout * omega_p(lambda, kappa, theta_k)^2 * (kappa^2 - omega_m(lambda, kappa, theta_k)^4) / (kappa^2 * omega_m(lambda, kappa, theta_k) * (omega_m(lambda, kappa, theta_k)^2 - omega_p(lambda, kappa, theta_k)^2) * besselj(1, rout * omega_m(lambda, kappa, theta_k))));
B1 = @(R, kappa, lambda, theta_k, rout) (8 * R^2 * (kappa^2 - omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2) * sqrt((kappa^2 - omega_m(lambda, kappa, theta_k)^4) * (kappa^2 - omega_p(lambda, kappa, theta_k)^4))) / ((rout^2 - 4 * R^2)^2 * (kappa^2 + omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2) * ...
(kappa^2 + omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2));
dr = @(r, lambda,kappa) ...
(2 * besselj(1, r(:,1) * omega_p(lambda, kappa, theta_k)) ./ ...
((omega_m(lambda, kappa, theta_k)^2 - ...
omega_p(lambda, kappa, theta_k)^2) * ...
(kappa^2 + omega_m(lambda, kappa, theta_k)^2 * ...
omega_p(lambda, kappa, theta_k)^2)) .* ...
((omega_m(lambda, kappa, theta_k)^2 * ...
(kappa^2 - omega_p(lambda, kappa, theta_k)^4)) ./ ...
(omega_p(lambda, kappa, theta_k) .* ...
(A1(R, kappa, lambda, theta_k, rout) + ...
(16 * R^2 * (kappa^2 - omega_m(lambda, kappa, theta_k)^2 * ...
omega_p(lambda, kappa, theta_k)^2))) ./ ...
(kappa^2 * omega_p(lambda, kappa, theta_k)^2 * ...
(rout^2 - 4 * R^2)^2) + (B1(R, kappa, lambda, theta_k, rout) ...
* omega_m(lambda, kappa, theta_k)^2 * ...
omega_p(lambda, kappa, theta_k) .* ...
sqrt((kappa^2 - omega_m(lambda, kappa, theta_k)^4) * ...
(kappa^2 - omega_p(lambda, kappa, theta_k)^4))) ./ kappa))) ...
- (2 * besselj(1, r(:,1) * omega_m(lambda, kappa, theta_k)) ./ ...
((omega_m(lambda, kappa, theta_k)^2 - ...
omega_p(lambda, kappa, theta_k)^2) * ...
(kappa^2 + omega_m(lambda, kappa, theta_k)^2 * ...
omega_p(lambda, kappa, theta_k)^2)) .* ...
((omega_p(lambda, kappa, theta_k)^2 * ...
(kappa^2 - omega_m(lambda, kappa, theta_k)^4)) ./ ...
(omega_m(lambda, kappa, theta_k) .* ...
(A1(R, kappa, lambda, theta_k, rout) + ...
(16 * R^2 * (kappa^2 - omega_m(lambda, kappa, theta_k)^2 * ...
omega_p(lambda, kappa, theta_k)^2)) ./ ...
(kappa^2 * omega_m(lambda, kappa, theta_k)^2 * ...
(rout^2 - 4 * R^2)^2) + (B1(R, kappa, lambda, theta_k, rout) ...
* omega_m(lambda, kappa, theta_k) * ...
omega_p(lambda, kappa, theta_k)^2 .* ...
sqrt((kappa^2 - omega_m(lambda, kappa, theta_k)^4) * ...
(kappa^2 - omega_p(lambda, kappa, theta_k)^4))) ./ kappa))) ...
- (16 * r(:,1) * R^2 * (omega_m(lambda, kappa, theta_k)^2 + ...
omega_p(lambda, kappa, theta_k)^2) .* ...
(kappa^2 - omega_m(lambda, kappa, theta_k)^2 * ...
omega_p(lambda, kappa, theta_k)^2)) ./ ...
(omega_p(lambda, kappa, theta_k)^2 * ...
omega_m(lambda, kappa, theta_k)^2 * ...
(rout^2 - 4 * R^2)^2 * (kappa^2 + ...
omega_m(lambda, kappa, theta_k)^2 * ...
omega_p(lambda, kappa, theta_k)^2)));
dtheta = @(r, lambda,kappa) ...
(2 * besselj(1, r(:,2) * omega_p(lambda, kappa, theta_k)) ./ ...
((kappa^2 + omega_m(lambda, kappa, theta_k)^2 * ...
omega_p(lambda, kappa, theta_k)^2) * ...
(omega_m(lambda, kappa, theta_k)^2 - ...
omega_p(lambda, kappa, theta_k)^2)) .* ...
(-sqrt((kappa^2 - omega_m(lambda, kappa, theta_k)^4) * ...
(kappa^2 - omega_p(lambda, kappa, theta_k)^4)) ./ ...
(omega_p(lambda, kappa, theta_k) .* ...
(A1(R, kappa, lambda, theta_k, rout) * kappa + ...
(16 * R^2 * (kappa^2 - omega_m(lambda, kappa, theta_k)^2 ...
* omega_p(lambda, kappa, theta_k)^2)) ./ ...
(kappa * omega_p(lambda, kappa, theta_k)^2 * ...
(rout^2 - 4 * R^2)^2) - B1(R, kappa, lambda, theta_k, rout) ...
* omega_p(lambda, kappa, theta_k) * ...
(kappa^2 - omega_m(lambda, kappa, theta_k)^4)))) ...
+ (2 * besselj(1, r(:,2) * omega_m(lambda, kappa, theta_k)) ./ ...
((kappa^2 + omega_m(lambda, kappa, theta_k)^2 * ...
omega_p(lambda, kappa, theta_k)^2) * ...
(omega_m(lambda, kappa, theta_k)^2 - ...
omega_p(lambda, kappa, theta_k)^2))) .* ...
(sqrt((kappa^2 - omega_m(lambda, kappa, theta_k)^4) * ...
(kappa^2 - omega_p(lambda, kappa, theta_k)^4)) ./ ...
(omega_m(lambda, kappa, theta_k) .* ...
(A1(R, kappa, lambda, theta_k, rout) * kappa + ...
(16 * R^2 * (kappa^2 - omega_m(lambda, kappa, theta_k)^2 * ...
omega_p(lambda, kappa, theta_k)^2)) ./ ...
(kappa * omega_m(lambda, kappa, theta_k)^2 * ...
(rout^2 - 4 * R^2)^2) + B1(R, kappa, lambda, theta_k, rout) * ...
omega_m(lambda, kappa, theta_k) * ...
(kappa^2 - omega_p(lambda, kappa, theta_k)^4)))) ...
+ (16 * r(:,2) * R^2 * (kappa^2 - omega_m(lambda, kappa, theta_k)^2 ...
* omega_p(lambda, kappa, theta_k)^2) * ...
sqrt((kappa^2 - omega_m(lambda, kappa, theta_k)^4) * ...
(kappa^2 - omega_p(lambda, kappa, theta_k)^4))) ./ ...
(kappa * omega_m(lambda, kappa, theta_k)^2 * ...
omega_p(lambda, kappa, theta_k)^2 * (rout^2 - 4 * R^2)^2 * ...
(kappa^2 + omega_m(lambda, kappa, theta_k)^2 * ...
omega_p(lambda, kappa, theta_k)^2)));
% Specify the maximum number of iterations for fitting
max_iterations = 1000000;
max_fun_evals = max_iterations;
% Fit the data for dr using lsqcurvefit
%options = optimoptions('lsqcurvefit', 'Algorithm', 'trust-region-reflective', 'MaxIterations', max_iterations,'MaxFunEvals',max_fun_evals);
%init = [dr([dr_data(:, 1),dtheta_data(:,1)],lambda_init,kappa_init),dtheta([dr_data(:, 1),dtheta_data(:,1)],lambda_init,kappa_init)];
options = optimoptions(@lsqcurvefit,'MaxIterations', max_iterations,'MaxFunEvals',max_fun_evals);
params = lsqcurvefit(@(params, r) [dr(r, params(1),params(2)),dtheta(r,params(1),params(2))], [lambda_init, kappa_init], [dr_data(:, 1),dtheta_data(:,1)], [dr_data(:, 2),dtheta_data(:,2)], [], [], options);
Local minimum possible.
lsqcurvefit stopped because the size of the current step is less than
the value of the step size tolerance.
% Extract the optimized parameters
lambda_fit = params(1)
lambda_fit = 2.2382
kappa_fit = params(2)
kappa_fit = 0.5899
% Use the initial guess for theta_k and R
theta_k_fit = theta_k;
R_fit = R;
% Plotting data
r_values = linspace(min(dr_data(:, 1)), max(dr_data(:, 1)), 1000).';
dr_fit = dr([r_values,r_values], params(1),params(2));
dtheta_fit = dtheta([r_values,r_values], params(1),params(2));
% Plotting
figure;
subplot(2, 1, 1);
plot(dr_data(:, 1), dr_data(:, 2), 'bo', 'DisplayName', 'Data');
hold on;
plot(r_values, dr_fit, 'r-', 'DisplayName', 'Fitted');
xlabel('r');
ylabel('dr');
title('Fitting dr');
legend;
subplot(2, 1, 2);
plot(dtheta_data(:, 1), dtheta_data(:, 2), 'bo', 'DisplayName', 'Data');
hold on;
plot(r_values, dtheta_fit, 'r-', 'DisplayName', 'Fitted');
xlabel('r');
ylabel('dtheta');
title('Fitting dtheta');
legend;
tuhin
am 2 Apr. 2024
⋮
-
Verknüpfen
Direkter Link zu diesem Kommentar
Bearbeitet: tuhin
am 2 Apr. 2024
Thanks, Torsten! I am facing one more problem. Can you please take a look and help me to fix it.
% 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];
% Manually specified initial guess values for parameters
lambda_init = 0.5; % Example initial guess for lambda
kappa_init = 0.5; % Example initial guess for kappa
theta_k_init = 0.5; % Example initial guess for theta_k
R_init = 250; % Example initial guess for R
% Constants
rout = 75; % Define rout
% Calculate omega_m and omega_p
omega_m = @(lambda, kappa, theta_k) sqrt(kappa / (2 * (lambda + 1))) * sqrt((lambda + 2) * cos(theta_k) - sqrt((lambda + 2)^2 * cos(theta_k)^2 - 4 * (lambda + 1)));
omega_p = @(lambda, kappa, theta_k) sqrt(kappa / (2 * (lambda + 1))) * sqrt((lambda + 2) * cos(theta_k) + sqrt((lambda + 2)^2 * cos(theta_k)^2 - 4 * (lambda + 1)));
% Define A1 and B1
A1 = @(R, kappa, lambda, theta_k, rout) (8 * R^2 * (kappa^2 - omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2)) / ((rout^2 - 4 * R^2)^2 * (kappa^2 + omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2)) * ...
(-2 * (omega_m(lambda, kappa, theta_k)^2 + omega_p(lambda, kappa, theta_k)^2) / (omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2) ...
+ rout * omega_m(lambda, kappa, theta_k)^2 * (kappa^2 - omega_p(lambda, kappa, theta_k)^4) / (kappa^2 * omega_p(lambda, kappa, theta_k) * (omega_m(lambda, kappa, theta_k)^2 - omega_p(lambda, kappa, theta_k)^2) * besselj(1, rout * omega_p(lambda, kappa, theta_k))) ...
- rout * omega_p(lambda, kappa, theta_k)^2 * (kappa^2 - omega_m(lambda, kappa, theta_k)^4) / (kappa^2 * omega_m(lambda, kappa, theta_k) * (omega_m(lambda, kappa, theta_k)^2 - omega_p(lambda, kappa, theta_k)^2) * besselj(1, rout * omega_m(lambda, kappa, theta_k))));
B1 = @(R, kappa, lambda, theta_k, rout)(8 * R^2 * (kappa^2 - omega_m(kappa, theta_k, lambda)^2 * omega_p(kappa, theta_k, lambda)^2) * sqrt((kappa^2 - omega_m(kappa, theta_k, lambda)^4) * (kappa^2 - omega_p(kappa, theta_k, lambda)^4))) / ...
((rout^2 - 4 * R^2)^2 * (kappa^2 + omega_m(kappa, theta_k, lambda)^2 * omega_p(kappa, theta_k, lambda)^2)) * ...
(2 / (omega_m(kappa, theta_k, lambda)^2 * omega_p(kappa, theta_k, lambda)^2 * kappa) + rout / (kappa * omega_m(kappa, theta_k, lambda) * (omega_m(kappa, theta_k, lambda)^2 - omega_p(kappa, theta_k, lambda)^2) * besselj(1, rout * omega_m(kappa, theta_k, lambda))) ...
- rout / (kappa * omega_p(kappa, theta_k, lambda) * (omega_m(kappa, theta_k, lambda)^2 - omega_p(kappa, theta_k, lambda)^2) * besselj(1, rout * omega_p(kappa, theta_k, lambda))));
dr = @(r, lambda,kappa,theta_k,R) ...
(2 * besselj(1, r(:,1) * omega_p(lambda, kappa, theta_k)) ./ ...
((omega_m(lambda, kappa, theta_k)^2 - ...
omega_p(lambda, kappa, theta_k)^2) * ...
(kappa^2 + omega_m(lambda, kappa, theta_k)^2 * ...
omega_p(lambda, kappa, theta_k)^2)) .* ...
((omega_m(lambda, kappa, theta_k)^2 * ...
(kappa^2 - omega_p(lambda, kappa, theta_k)^4)) ./ ...
(omega_p(lambda, kappa, theta_k) .* ...
(A1(R, kappa, lambda, theta_k, rout) + ...
(16 * R^2 * (kappa^2 - omega_m(lambda, kappa, theta_k)^2 * ...
omega_p(lambda, kappa, theta_k)^2))) ./ ...
(kappa^2 * omega_p(lambda, kappa, theta_k)^2 * ...
(rout^2 - 4 * R^2)^2) + (B1(R, kappa, lambda, theta_k, rout) ...
* omega_m(lambda, kappa, theta_k)^2 * ...
omega_p(lambda, kappa, theta_k) .* ...
sqrt((kappa^2 - omega_m(lambda, kappa, theta_k)^4) * ...
(kappa^2 - omega_p(lambda, kappa, theta_k)^4))) ./ kappa))) ...
- (2 * besselj(1, r(:,1) * omega_m(lambda, kappa, theta_k)) ./ ...
((omega_m(lambda, kappa, theta_k)^2 - ...
omega_p(lambda, kappa, theta_k)^2) * ...
(kappa^2 + omega_m(lambda, kappa, theta_k)^2 * ...
omega_p(lambda, kappa, theta_k)^2)) .* ...
((omega_p(lambda, kappa, theta_k)^2 * ...
(kappa^2 - omega_m(lambda, kappa, theta_k)^4)) ./ ...
(omega_m(lambda, kappa, theta_k) .* ...
(A1(R, kappa, lambda, theta_k, rout) + ...
(16 * R^2 * (kappa^2 - omega_m(lambda, kappa, theta_k)^2 * ...
omega_p(lambda, kappa, theta_k)^2)) ./ ...
(kappa^2 * omega_m(lambda, kappa, theta_k)^2 * ...
(rout^2 - 4 * R^2)^2) + (B1(R, kappa, lambda, theta_k, rout) ...
* omega_m(lambda, kappa, theta_k) * ...
omega_p(lambda, kappa, theta_k)^2 .* ...
sqrt((kappa^2 - omega_m(lambda, kappa, theta_k)^4) * ...
(kappa^2 - omega_p(lambda, kappa, theta_k)^4))) ./ kappa))) ...
- (16 * r(:,1) * R^2 * (omega_m(lambda, kappa, theta_k)^2 + ...
omega_p(lambda, kappa, theta_k)^2) .* ...
(kappa^2 - omega_m(lambda, kappa, theta_k)^2 * ...
omega_p(lambda, kappa, theta_k)^2)) ./ ...
(omega_p(lambda, kappa, theta_k)^2 * ...
omega_m(lambda, kappa, theta_k)^2 * ...
(rout^2 - 4 * R^2)^2 * (kappa^2 + ...
omega_m(lambda, kappa, theta_k)^2 * ...
omega_p(lambda, kappa, theta_k)^2)));
dtheta = @(r, lambda,kappa,theta_k,R) ...
(2 * besselj(1, r(:,2) * omega_p(lambda, kappa, theta_k)) ./ ...
((kappa^2 + omega_m(lambda, kappa, theta_k)^2 * ...
omega_p(lambda, kappa, theta_k)^2) * ...
(omega_m(lambda, kappa, theta_k)^2 - ...
omega_p(lambda, kappa, theta_k)^2)) .* ...
(-sqrt((kappa^2 - omega_m(lambda, kappa, theta_k)^4) * ...
(kappa^2 - omega_p(lambda, kappa, theta_k)^4)) ./ ...
(omega_p(lambda, kappa, theta_k) .* ...
(A1(R, kappa, lambda, theta_k, rout) * kappa + ...
(16 * R^2 * (kappa^2 - omega_m(lambda, kappa, theta_k)^2 ...
* omega_p(lambda, kappa, theta_k)^2)) ./ ...
(kappa * omega_p(lambda, kappa, theta_k)^2 * ...
(rout^2 - 4 * R^2)^2) - B1(R, kappa, lambda, theta_k, rout) ...
* omega_p(lambda, kappa, theta_k) * ...
(kappa^2 - omega_m(lambda, kappa, theta_k)^4)))) ...
+ (2 * besselj(1, r(:,2) * omega_m(lambda, kappa, theta_k)) ./ ...
((kappa^2 + omega_m(lambda, kappa, theta_k)^2 * ...
omega_p(lambda, kappa, theta_k)^2) * ...
(omega_m(lambda, kappa, theta_k)^2 - ...
omega_p(lambda, kappa, theta_k)^2))) .* ...
(sqrt((kappa^2 - omega_m(lambda, kappa, theta_k)^4) * ...
(kappa^2 - omega_p(lambda, kappa, theta_k)^4)) ./ ...
(omega_m(lambda, kappa, theta_k) .* ...
(A1(R, kappa, lambda, theta_k, rout) * kappa + ...
(16 * R^2 * (kappa^2 - omega_m(lambda, kappa, theta_k)^2 * ...
omega_p(lambda, kappa, theta_k)^2)) ./ ...
(kappa * omega_m(lambda, kappa, theta_k)^2 * ...
(rout^2 - 4 * R^2)^2) + B1(R, kappa, lambda, theta_k, rout) * ...
omega_m(lambda, kappa, theta_k) * ...
(kappa^2 - omega_p(lambda, kappa, theta_k)^4)))) ...
+ (16 * r(:,2) * R^2 * (kappa^2 - omega_m(lambda, kappa, theta_k)^2 ...
* omega_p(lambda, kappa, theta_k)^2) * ...
sqrt((kappa^2 - omega_m(lambda, kappa, theta_k)^4) * ...
(kappa^2 - omega_p(lambda, kappa, theta_k)^4))) ./ ...
(kappa * omega_m(lambda, kappa, theta_k)^2 * ...
omega_p(lambda, kappa, theta_k)^2 * (rout^2 - 4 * R^2)^2 * ...
(kappa^2 + omega_m(lambda, kappa, theta_k)^2 * ...
omega_p(lambda, kappa, theta_k)^2)));
% Specify the maximum number of iterations for fitting
max_iterations = 1000000;
max_fun_evals = max_iterations;
% Fit the data for dr using lsqcurvefit
%options = optimoptions('lsqcurvefit', 'Algorithm', 'trust-region-reflective', 'MaxIterations', max_iterations,'MaxFunEvals',max_fun_evals);
%init = [dr([dr_data(:, 1),dtheta_data(:,1)],lambda_init,kappa_init),dtheta([dr_data(:, 1),dtheta_data(:,1)],lambda_init,kappa_init)];
options = optimoptions(@lsqcurvefit,'MaxIterations', max_iterations,'MaxFunEvals',max_fun_evals);
params = lsqcurvefit(@(params, r) [dr(r, params(1),params(2),params(3),params(4)),dtheta(r,params(1),params(2),params(3),params(4))], [lambda_init, kappa_init,theta_k_init,R_init], [dr_data(:, 1),dtheta_data(:,1)], [dr_data(:, 2),dtheta_data(:,2)], [], [], options);
% Extract the optimized parameters
lambda_fit = params(1)
kappa_fit = params(2)
theta_k_fit = params(3)
R_fit = params(4)
% Use the initial guess for theta_k and R
%theta_k_fit = theta_k;
%R_fit = R;
% Plotting data
r_values = linspace(min(dr_data(:, 1)), max(dr_data(:, 1)), 1000).';
dr_fit = dr([r_values,r_values,r_values,r_values], params(1),params(2),params(3),params(4));
dtheta_fit = dtheta([r_values,r_values,r_values,r_values], params(1),params(2),params(3),params(4));
% Plotting
figure;
subplot(2, 1, 1);
plot(dr_data(:, 1), dr_data(:, 2), 'bo', 'DisplayName', 'Data');
hold on;
plot(r_values, dr_fit, 'r-', 'DisplayName', 'Fitted');
xlabel('r');
ylabel('dr');
title('Fitting dr');
legend;
subplot(2, 1, 2);
plot(dtheta_data(:, 1), dtheta_data(:, 2), 'bo', 'DisplayName', 'Data');
hold on;
plot(r_values, dtheta_fit, 'r-', 'DisplayName', 'Fitted');
xlabel('r');
ylabel('dtheta');
title('Fitting dtheta');
legend;
Here I am trying to fit the data using 4 parameters and getting some imaginary part as the fitted values. Please check those values:
lambda_fit =
0.5683 - 0.0093i
kappa_fit =
0.4829 - 0.0029i
theta_k_fit =
-1.1943e-05 + 5.5348e-04i
R_fit =
18.5555 + 0.1024i
Please let me know how to do it.
One more thing I also want the R_fit values should be anythong above 100.
One more thing I also want the theta_k_fit values should be in between 0 to 1.57.
Please take a look and let me know how should i modify this code.
tuhin
am 2 Apr. 2024
⋮
-
Verknüpfen
Direkter Link zu diesem Kommentar
Bearbeitet: Torsten
am 2 Apr. 2024
I tried to solve this with this modification but it is not even working. Please see below:
% 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];
% Manually specified initial guess values for parameters
lambda_init = 0.5; % Example initial guess for lambda
kappa_init = 0.5; % Example initial guess for kappa
theta_k_init = 0.5; % Example initial guess for theta_k
R_init = 150; % Example initial guess for R
% Constants
rout = 75; % Define rout
% Calculate omega_m and omega_p
omega_m = @(lambda, kappa, theta_k) sqrt(kappa / (2 * (lambda + 1))) * sqrt((lambda + 2) * cos(theta_k) - sqrt((lambda + 2)^2 * cos(theta_k)^2 - 4 * (lambda + 1)));
omega_p = @(lambda, kappa, theta_k) sqrt(kappa / (2 * (lambda + 1))) * sqrt((lambda + 2) * cos(theta_k) + sqrt((lambda + 2)^2 * cos(theta_k)^2 - 4 * (lambda + 1)));
% Define A1 and B1
A1 = @(R, kappa, lambda, theta_k, rout) (8 * R^2 * (kappa^2 - omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2)) / ((rout^2 - 4 * R^2)^2 * (kappa^2 + omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2)) * ...
(-2 * (omega_m(lambda, kappa, theta_k)^2 + omega_p(lambda, kappa, theta_k)^2) / (omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2) ...
+ rout * omega_m(lambda, kappa, theta_k)^2 * (kappa^2 - omega_p(lambda, kappa, theta_k)^4) / (kappa^2 * omega_p(lambda, kappa, theta_k) * (omega_m(lambda, kappa, theta_k)^2 - omega_p(lambda, kappa, theta_k)^2) * besselj(1, rout * omega_p(lambda, kappa, theta_k))) ...
- rout * omega_p(lambda, kappa, theta_k)^2 * (kappa^2 - omega_m(lambda, kappa, theta_k)^4) / (kappa^2 * omega_m(lambda, kappa, theta_k) * (omega_m(lambda, kappa, theta_k)^2 - omega_p(lambda, kappa, theta_k)^2) * besselj(1, rout * omega_m(lambda, kappa, theta_k))));
B1 = @(R, kappa, lambda, theta_k, rout)(8 * R^2 * (kappa^2 - omega_m(kappa, theta_k, lambda)^2 * omega_p(kappa, theta_k, lambda)^2) * sqrt((kappa^2 - omega_m(kappa, theta_k, lambda)^4) * (kappa^2 - omega_p(kappa, theta_k, lambda)^4))) / ...
((rout^2 - 4 * R^2)^2 * (kappa^2 + omega_m(kappa, theta_k, lambda)^2 * omega_p(kappa, theta_k, lambda)^2)) * ...
(2 / (omega_m(kappa, theta_k, lambda)^2 * omega_p(kappa, theta_k, lambda)^2 * kappa) + rout / (kappa * omega_m(kappa, theta_k, lambda) * (omega_m(kappa, theta_k, lambda)^2 - omega_p(kappa, theta_k, lambda)^2) * besselj(1, rout * omega_m(kappa, theta_k, lambda))) ...
- rout / (kappa * omega_p(kappa, theta_k, lambda) * (omega_m(kappa, theta_k, lambda)^2 - omega_p(kappa, theta_k, lambda)^2) * besselj(1, rout * omega_p(kappa, theta_k, lambda))));
dr = @(r, lambda,kappa,theta_k,R) ...
(2 * besselj(1, r(:,1) * omega_p(lambda, kappa, theta_k)) ./ ...
((omega_m(lambda, kappa, theta_k)^2 - ...
omega_p(lambda, kappa, theta_k)^2) * ...
(kappa^2 + omega_m(lambda, kappa, theta_k)^2 * ...
omega_p(lambda, kappa, theta_k)^2)) .* ...
((omega_m(lambda, kappa, theta_k)^2 * ...
(kappa^2 - omega_p(lambda, kappa, theta_k)^4)) ./ ...
(omega_p(lambda, kappa, theta_k) .* ...
(A1(R, kappa, lambda, theta_k, rout) + ...
(16 * R^2 * (kappa^2 - omega_m(lambda, kappa, theta_k)^2 * ...
omega_p(lambda, kappa, theta_k)^2))) ./ ...
(kappa^2 * omega_p(lambda, kappa, theta_k)^2 * ...
(rout^2 - 4 * R^2)^2) + (B1(R, kappa, lambda, theta_k, rout) ...
* omega_m(lambda, kappa, theta_k)^2 * ...
omega_p(lambda, kappa, theta_k) .* ...
sqrt((kappa^2 - omega_m(lambda, kappa, theta_k)^4) * ...
(kappa^2 - omega_p(lambda, kappa, theta_k)^4))) ./ kappa))) ...
- (2 * besselj(1, r(:,1) * omega_m(lambda, kappa, theta_k)) ./ ...
((omega_m(lambda, kappa, theta_k)^2 - ...
omega_p(lambda, kappa, theta_k)^2) * ...
(kappa^2 + omega_m(lambda, kappa, theta_k)^2 * ...
omega_p(lambda, kappa, theta_k)^2)) .* ...
((omega_p(lambda, kappa, theta_k)^2 * ...
(kappa^2 - omega_m(lambda, kappa, theta_k)^4)) ./ ...
(omega_m(lambda, kappa, theta_k) .* ...
(A1(R, kappa, lambda, theta_k, rout) + ...
(16 * R^2 * (kappa^2 - omega_m(lambda, kappa, theta_k)^2 * ...
omega_p(lambda, kappa, theta_k)^2)) ./ ...
(kappa^2 * omega_m(lambda, kappa, theta_k)^2 * ...
(rout^2 - 4 * R^2)^2) + (B1(R, kappa, lambda, theta_k, rout) ...
* omega_m(lambda, kappa, theta_k) * ...
omega_p(lambda, kappa, theta_k)^2 .* ...
sqrt((kappa^2 - omega_m(lambda, kappa, theta_k)^4) * ...
(kappa^2 - omega_p(lambda, kappa, theta_k)^4))) ./ kappa))) ...
- (16 * r(:,1) * R^2 * (omega_m(lambda, kappa, theta_k)^2 + ...
omega_p(lambda, kappa, theta_k)^2) .* ...
(kappa^2 - omega_m(lambda, kappa, theta_k)^2 * ...
omega_p(lambda, kappa, theta_k)^2)) ./ ...
(omega_p(lambda, kappa, theta_k)^2 * ...
omega_m(lambda, kappa, theta_k)^2 * ...
(rout^2 - 4 * R^2)^2 * (kappa^2 + ...
omega_m(lambda, kappa, theta_k)^2 * ...
omega_p(lambda, kappa, theta_k)^2)));
dtheta = @(r, lambda,kappa,theta_k,R) ...
(2 * besselj(1, r(:,2) * omega_p(lambda, kappa, theta_k)) ./ ...
((kappa^2 + omega_m(lambda, kappa, theta_k)^2 * ...
omega_p(lambda, kappa, theta_k)^2) * ...
(omega_m(lambda, kappa, theta_k)^2 - ...
omega_p(lambda, kappa, theta_k)^2)) .* ...
(-sqrt((kappa^2 - omega_m(lambda, kappa, theta_k)^4) * ...
(kappa^2 - omega_p(lambda, kappa, theta_k)^4)) ./ ...
(omega_p(lambda, kappa, theta_k) .* ...
(A1(R, kappa, lambda, theta_k, rout) * kappa + ...
(16 * R^2 * (kappa^2 - omega_m(lambda, kappa, theta_k)^2 ...
* omega_p(lambda, kappa, theta_k)^2)) ./ ...
(kappa * omega_p(lambda, kappa, theta_k)^2 * ...
(rout^2 - 4 * R^2)^2) - B1(R, kappa, lambda, theta_k, rout) ...
* omega_p(lambda, kappa, theta_k) * ...
(kappa^2 - omega_m(lambda, kappa, theta_k)^4)))) ...
+ (2 * besselj(1, r(:,2) * omega_m(lambda, kappa, theta_k)) ./ ...
((kappa^2 + omega_m(lambda, kappa, theta_k)^2 * ...
omega_p(lambda, kappa, theta_k)^2) * ...
(omega_m(lambda, kappa, theta_k)^2 - ...
omega_p(lambda, kappa, theta_k)^2))) .* ...
(sqrt((kappa^2 - omega_m(lambda, kappa, theta_k)^4) * ...
(kappa^2 - omega_p(lambda, kappa, theta_k)^4)) ./ ...
(omega_m(lambda, kappa, theta_k) .* ...
(A1(R, kappa, lambda, theta_k, rout) * kappa + ...
(16 * R^2 * (kappa^2 - omega_m(lambda, kappa, theta_k)^2 * ...
omega_p(lambda, kappa, theta_k)^2)) ./ ...
(kappa * omega_m(lambda, kappa, theta_k)^2 * ...
(rout^2 - 4 * R^2)^2) + B1(R, kappa, lambda, theta_k, rout) * ...
omega_m(lambda, kappa, theta_k) * ...
(kappa^2 - omega_p(lambda, kappa, theta_k)^4)))) ...
+ (16 * r(:,2) * R^2 * (kappa^2 - omega_m(lambda, kappa, theta_k)^2 ...
* omega_p(lambda, kappa, theta_k)^2) * ...
sqrt((kappa^2 - omega_m(lambda, kappa, theta_k)^4) * ...
(kappa^2 - omega_p(lambda, kappa, theta_k)^4))) ./ ...
(kappa * omega_m(lambda, kappa, theta_k)^2 * ...
omega_p(lambda, kappa, theta_k)^2 * (rout^2 - 4 * R^2)^2 * ...
(kappa^2 + omega_m(lambda, kappa, theta_k)^2 * ...
omega_p(lambda, kappa, theta_k)^2)));
% Specify the maximum number of iterations for fitting
max_iterations = 1000000;
max_fun_evals = max_iterations;
% Fit the data for dr using lsqcurvefit
options = optimoptions(@lsqcurvefit,'MaxIterations', max_iterations,'MaxFunctionEvaluations',max_fun_evals);
lb = [0, 0, 0, 100]; % Lower bounds for lambda, kappa, theta_k, and R
ub = [Inf, Inf, 1.57, Inf]; % Upper bounds for lambda, kappa, theta_k, and R
init = [dr([dr_data(:, 1),dtheta_data(:,1)],lambda_init,kappa_init,theta_k_init,R_init),dtheta([dr_data(:, 1),dtheta_data(:,1)],lambda_init,kappa_init,theta_k_init,R_init)]
init =
1.0e+09 *
-0.0000 - 0.0000i 0.0005 + 0.0000i
0.0000 - 0.0000i 0.0005 + 0.0001i
0.0000 - 0.0000i -0.0001 + 0.0001i
-0.0000 + 0.0000i -0.0006 - 0.0000i
-0.0001 + 0.0001i -0.0002 - 0.0002i
-0.0000 - 0.0000i 0.0008 - 0.0001i
0.0001 - 0.0001i 0.0007 + 0.0003i
0.0001 - 0.0000i -0.0008 + 0.0003i
-0.0001 + 0.0002i -0.0016 - 0.0004i
-0.0002 + 0.0001i 0.0005 - 0.0007i
0.0000 - 0.0002i 0.0031 + 0.0002i
0.0004 - 0.0003i 0.0009 + 0.0013i
0.0002 + 0.0002i -0.0047 + 0.0004i
-0.0006 + 0.0007i -0.0044 - 0.0020i
-0.0008 + 0.0001i 0.0058 - 0.0018i
0.0007 - 0.0011i 0.0111 + 0.0025i
0.0017 - 0.0008i -0.0037 + 0.0046i
-0.0002 + 0.0015i -0.0219 - 0.0017i
-0.0032 + 0.0022i -0.0062 - 0.0092i
-0.0016 - 0.0014i 0.0356 - 0.0023i
0.0048 - 0.0049i 0.0320 + 0.0152i
0.0057 - 0.0003i -0.0452 + 0.0131i
-0.0053 + 0.0086i -0.0842 - 0.0196i
-0.0133 + 0.0056i 0.0320 - 0.0353i
0.0018 - 0.0123i 0.1711 + 0.0148i
0.0252 - 0.0173i 0.0428 + 0.0726i
0.0116 + 0.0120i -0.2852 + 0.0160i
-0.0390 + 0.0386i -0.2459 - 0.1224i
-0.0441 + 0.0009i 0.3736 - 0.1013i
0.0450 - 0.0697i 0.6682 + 0.1631i
0.1069 - 0.0427i -0.2866 + 0.2809i
-0.0190 + 0.1019i -1.3887 - 0.1317i
-0.2065 + 0.1383i -0.2970 - 0.5910i
-0.0874 - 0.1046i 2.3639 - 0.1084i
params = lsqcurvefit(@(params, r) [dr(r, params(1),params(2),params(3),params(4)),dtheta(r,params(1),params(2),params(3),params(4))], [lambda_init, kappa_init,theta_k_init,R_init], [dr_data(:, 1),dtheta_data(:,1)], [dr_data(:, 2),dtheta_data(:,2)], lb, ub, options);
Error using lsqncommon (line 65)
Lower and upper bounds not supported with complex-valued initial function or Jacobian evaluation.
Lower and upper bounds not supported with complex-valued initial function or Jacobian evaluation.
Error in lsqcurvefit (line 306)
lsqncommon(funfcn,xCurrent,lb,ub,options,defaultopt,optimgetFlag,caller,...
% Extract the optimized parameters
lambda_fit = params(1);
kappa_fit = params(2);
theta_k_fit = params(3);
R_fit = params(4);
% Use the initial guess for theta_k and R
%theta_k_fit = theta_k;
%R_fit = R;
% Plotting data
r_values = linspace(min(dr_data(:, 1)), max(dr_data(:, 1)), 1000).';
dr_fit = dr([r_values,r_values], params(1),params(2),params(3),params(4));
dtheta_fit = dtheta([r_values,r_values], params(1),params(2),params(3),params(4));
% Plotting
figure;
subplot(2, 1, 1);
plot(dr_data(:, 1), dr_data(:, 2), 'bo', 'DisplayName', 'Data');
hold on;
plot(r_values, dr_fit, 'r-', 'DisplayName', 'Fitted');
xlabel('r');
ylabel('dr');
title('Fitting dr');
legend;
subplot(2, 1, 2);
plot(dtheta_data(:, 1), dtheta_data(:, 2), 'bo', 'DisplayName', 'Data');
hold on;
plot(r_values, dtheta_fit, 'r-', 'DisplayName', 'Fitted');
xlabel('r');
ylabel('dtheta');
title('Fitting dtheta');
legend;
Torsten
am 2 Apr. 2024
⋮
-
Verknüpfen
Direkter Link zu diesem Kommentar
Bearbeitet: Torsten
am 2 Apr. 2024
Your function produces complex numbers (most probably because of a parameter constellation for which you take the squareroot of a negative number) (see above). But setting bounds on complex numbers (lb, ub) doesn't make sense.
Since this already happens for the initial values, you should adjust them such that the squareroots in your functions for dr and dtheta produce real values. But it's not guaranteed that during the course of the iteration, complex numbers could again be produced. You could try to avoid this be returning high values for the residuals in such a case.
tuhin
am 2 Apr. 2024
⋮
-
Verknüpfen
Direkter Link zu diesem Kommentar
I have no idea either. I would really appreciate it if someone could assist me with this.
Sam Chak
am 2 Apr. 2024
⋮
-
Verknüpfen
Direkter Link zu diesem Kommentar
Hi @tuhin, you can click on the Indentation icon 
to copy/paste your code, and then click on the Play icon 
to run the code on MATLAB Answers forum. This allows to you read the error messages if there are bugs in the code.
to copy/paste your code, and then click on the Play icon to run the code on MATLAB Answers forum. This allows to you read the error messages if there are bugs in the code.
tuhin
am 2 Apr. 2024
⋮
-
Verknüpfen
Direkter Link zu diesem Kommentar
% 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];
% Manually specified initial guess values for parameters
lambda_init = 0.5; % Example initial guess for lambda
kappa_init = 0.5; % Example initial guess for kappa
theta_k_init = 0.5; % Example initial guess for theta_k
R_init = 150; % Example initial guess for R
% Constants
rout = 75; % Define rout
% Calculate omega_m and omega_p
omega_m = @(lambda, kappa, theta_k) sqrt(kappa / (2 * (lambda + 1))) * sqrt((lambda + 2) * cos(theta_k) - sqrt((lambda + 2)^2 * cos(theta_k)^2 - 4 * (lambda + 1)));
omega_p = @(lambda, kappa, theta_k) sqrt(kappa / (2 * (lambda + 1))) * sqrt((lambda + 2) * cos(theta_k) + sqrt((lambda + 2)^2 * cos(theta_k)^2 - 4 * (lambda + 1)));
% Define A1 and B1
A1 = @(R, kappa, lambda, theta_k, rout) (8 * R^2 * (kappa^2 - omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2)) / ((rout^2 - 4 * R^2)^2 * (kappa^2 + omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2)) * ...
(-2 * (omega_m(lambda, kappa, theta_k)^2 + omega_p(lambda, kappa, theta_k)^2) / (omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2) ...
+ rout * omega_m(lambda, kappa, theta_k)^2 * (kappa^2 - omega_p(lambda, kappa, theta_k)^4) / (kappa^2 * omega_p(lambda, kappa, theta_k) * (omega_m(lambda, kappa, theta_k)^2 - omega_p(lambda, kappa, theta_k)^2) * besselj(1, rout * omega_p(lambda, kappa, theta_k))) ...
- rout * omega_p(lambda, kappa, theta_k)^2 * (kappa^2 - omega_m(lambda, kappa, theta_k)^4) / (kappa^2 * omega_m(lambda, kappa, theta_k) * (omega_m(lambda, kappa, theta_k)^2 - omega_p(lambda, kappa, theta_k)^2) * besselj(1, rout * omega_m(lambda, kappa, theta_k))));
B1 = @(R, kappa, lambda, theta_k, rout)(8 * R^2 * (kappa^2 - omega_m(kappa, theta_k, lambda)^2 * omega_p(kappa, theta_k, lambda)^2) * sqrt((kappa^2 - omega_m(kappa, theta_k, lambda)^4) * (kappa^2 - omega_p(kappa, theta_k, lambda)^4))) / ...
((rout^2 - 4 * R^2)^2 * (kappa^2 + omega_m(kappa, theta_k, lambda)^2 * omega_p(kappa, theta_k, lambda)^2)) * ...
(2 / (omega_m(kappa, theta_k, lambda)^2 * omega_p(kappa, theta_k, lambda)^2 * kappa) + rout / (kappa * omega_m(kappa, theta_k, lambda) * (omega_m(kappa, theta_k, lambda)^2 - omega_p(kappa, theta_k, lambda)^2) * besselj(1, rout * omega_m(kappa, theta_k, lambda))) ...
- rout / (kappa * omega_p(kappa, theta_k, lambda) * (omega_m(kappa, theta_k, lambda)^2 - omega_p(kappa, theta_k, lambda)^2) * besselj(1, rout * omega_p(kappa, theta_k, lambda))));
dr = @(r, lambda,kappa,theta_k,R) ...
(2 * besselj(1, r(:,1) * omega_p(lambda, kappa, theta_k)) ./ ...
((omega_m(lambda, kappa, theta_k)^2 - ...
omega_p(lambda, kappa, theta_k)^2) * ...
(kappa^2 + omega_m(lambda, kappa, theta_k)^2 * ...
omega_p(lambda, kappa, theta_k)^2)) .* ...
((omega_m(lambda, kappa, theta_k)^2 * ...
(kappa^2 - omega_p(lambda, kappa, theta_k)^4)) ./ ...
(omega_p(lambda, kappa, theta_k) .* ...
(A1(R, kappa, lambda, theta_k, rout) + ...
(16 * R^2 * (kappa^2 - omega_m(lambda, kappa, theta_k)^2 * ...
omega_p(lambda, kappa, theta_k)^2))) ./ ...
(kappa^2 * omega_p(lambda, kappa, theta_k)^2 * ...
(rout^2 - 4 * R^2)^2) + (B1(R, kappa, lambda, theta_k, rout) ...
* omega_m(lambda, kappa, theta_k)^2 * ...
omega_p(lambda, kappa, theta_k) .* ...
sqrt((kappa^2 - omega_m(lambda, kappa, theta_k)^4) * ...
(kappa^2 - omega_p(lambda, kappa, theta_k)^4))) ./ kappa))) ...
- (2 * besselj(1, r(:,1) * omega_m(lambda, kappa, theta_k)) ./ ...
((omega_m(lambda, kappa, theta_k)^2 - ...
omega_p(lambda, kappa, theta_k)^2) * ...
(kappa^2 + omega_m(lambda, kappa, theta_k)^2 * ...
omega_p(lambda, kappa, theta_k)^2)) .* ...
((omega_p(lambda, kappa, theta_k)^2 * ...
(kappa^2 - omega_m(lambda, kappa, theta_k)^4)) ./ ...
(omega_m(lambda, kappa, theta_k) .* ...
(A1(R, kappa, lambda, theta_k, rout) + ...
(16 * R^2 * (kappa^2 - omega_m(lambda, kappa, theta_k)^2 * ...
omega_p(lambda, kappa, theta_k)^2)) ./ ...
(kappa^2 * omega_m(lambda, kappa, theta_k)^2 * ...
(rout^2 - 4 * R^2)^2) + (B1(R, kappa, lambda, theta_k, rout) ...
* omega_m(lambda, kappa, theta_k) * ...
omega_p(lambda, kappa, theta_k)^2 .* ...
sqrt((kappa^2 - omega_m(lambda, kappa, theta_k)^4) * ...
(kappa^2 - omega_p(lambda, kappa, theta_k)^4))) ./ kappa))) ...
- (16 * r(:,1) * R^2 * (omega_m(lambda, kappa, theta_k)^2 + ...
omega_p(lambda, kappa, theta_k)^2) .* ...
(kappa^2 - omega_m(lambda, kappa, theta_k)^2 * ...
omega_p(lambda, kappa, theta_k)^2)) ./ ...
(omega_p(lambda, kappa, theta_k)^2 * ...
omega_m(lambda, kappa, theta_k)^2 * ...
(rout^2 - 4 * R^2)^2 * (kappa^2 + ...
omega_m(lambda, kappa, theta_k)^2 * ...
omega_p(lambda, kappa, theta_k)^2)));
dtheta = @(r, lambda,kappa,theta_k,R) ...
(2 * besselj(1, r(:,2) * omega_p(lambda, kappa, theta_k)) ./ ...
((kappa^2 + omega_m(lambda, kappa, theta_k)^2 * ...
omega_p(lambda, kappa, theta_k)^2) * ...
(omega_m(lambda, kappa, theta_k)^2 - ...
omega_p(lambda, kappa, theta_k)^2)) .* ...
(-sqrt((kappa^2 - omega_m(lambda, kappa, theta_k)^4) * ...
(kappa^2 - omega_p(lambda, kappa, theta_k)^4)) ./ ...
(omega_p(lambda, kappa, theta_k) .* ...
(A1(R, kappa, lambda, theta_k, rout) * kappa + ...
(16 * R^2 * (kappa^2 - omega_m(lambda, kappa, theta_k)^2 ...
* omega_p(lambda, kappa, theta_k)^2)) ./ ...
(kappa * omega_p(lambda, kappa, theta_k)^2 * ...
(rout^2 - 4 * R^2)^2) - B1(R, kappa, lambda, theta_k, rout) ...
* omega_p(lambda, kappa, theta_k) * ...
(kappa^2 - omega_m(lambda, kappa, theta_k)^4)))) ...
+ (2 * besselj(1, r(:,2) * omega_m(lambda, kappa, theta_k)) ./ ...
((kappa^2 + omega_m(lambda, kappa, theta_k)^2 * ...
omega_p(lambda, kappa, theta_k)^2) * ...
(omega_m(lambda, kappa, theta_k)^2 - ...
omega_p(lambda, kappa, theta_k)^2))) .* ...
(sqrt((kappa^2 - omega_m(lambda, kappa, theta_k)^4) * ...
(kappa^2 - omega_p(lambda, kappa, theta_k)^4)) ./ ...
(omega_m(lambda, kappa, theta_k) .* ...
(A1(R, kappa, lambda, theta_k, rout) * kappa + ...
(16 * R^2 * (kappa^2 - omega_m(lambda, kappa, theta_k)^2 * ...
omega_p(lambda, kappa, theta_k)^2)) ./ ...
(kappa * omega_m(lambda, kappa, theta_k)^2 * ...
(rout^2 - 4 * R^2)^2) + B1(R, kappa, lambda, theta_k, rout) * ...
omega_m(lambda, kappa, theta_k) * ...
(kappa^2 - omega_p(lambda, kappa, theta_k)^4)))) ...
+ (16 * r(:,2) * R^2 * (kappa^2 - omega_m(lambda, kappa, theta_k)^2 ...
* omega_p(lambda, kappa, theta_k)^2) * ...
sqrt((kappa^2 - omega_m(lambda, kappa, theta_k)^4) * ...
(kappa^2 - omega_p(lambda, kappa, theta_k)^4))) ./ ...
(kappa * omega_m(lambda, kappa, theta_k)^2 * ...
omega_p(lambda, kappa, theta_k)^2 * (rout^2 - 4 * R^2)^2 * ...
(kappa^2 + omega_m(lambda, kappa, theta_k)^2 * ...
omega_p(lambda, kappa, theta_k)^2)));
% Specify the maximum number of iterations for fitting
max_iterations = 1000000;
max_fun_evals = max_iterations;
% Fit the data for dr using lsqcurvefit
options = optimoptions(@lsqcurvefit,'MaxIterations', max_iterations,'MaxFunctionEvaluations',max_fun_evals);
lb = [0, 0, 0, 100]; % Lower bounds for lambda, kappa, theta_k, and R
ub = [Inf, Inf, 1.57, Inf]; % Upper bounds for lambda, kappa, theta_k, and R
init = [dr([dr_data(:, 1),dtheta_data(:,1)],lambda_init,kappa_init,theta_k_init,R_init),dtheta([dr_data(:, 1),dtheta_data(:,1)],lambda_init,kappa_init,theta_k_init,R_init)]
init =
1.0e+09 *
-0.0000 - 0.0000i 0.0005 + 0.0000i
0.0000 - 0.0000i 0.0005 + 0.0001i
0.0000 - 0.0000i -0.0001 + 0.0001i
-0.0000 + 0.0000i -0.0006 - 0.0000i
-0.0001 + 0.0001i -0.0002 - 0.0002i
-0.0000 - 0.0000i 0.0008 - 0.0001i
0.0001 - 0.0001i 0.0007 + 0.0003i
0.0001 - 0.0000i -0.0008 + 0.0003i
-0.0001 + 0.0002i -0.0016 - 0.0004i
-0.0002 + 0.0001i 0.0005 - 0.0007i
0.0000 - 0.0002i 0.0031 + 0.0002i
0.0004 - 0.0003i 0.0009 + 0.0013i
0.0002 + 0.0002i -0.0047 + 0.0004i
-0.0006 + 0.0007i -0.0044 - 0.0020i
-0.0008 + 0.0001i 0.0058 - 0.0018i
0.0007 - 0.0011i 0.0111 + 0.0025i
0.0017 - 0.0008i -0.0037 + 0.0046i
-0.0002 + 0.0015i -0.0219 - 0.0017i
-0.0032 + 0.0022i -0.0062 - 0.0092i
-0.0016 - 0.0014i 0.0356 - 0.0023i
0.0048 - 0.0049i 0.0320 + 0.0152i
0.0057 - 0.0003i -0.0452 + 0.0131i
-0.0053 + 0.0086i -0.0842 - 0.0196i
-0.0133 + 0.0056i 0.0320 - 0.0353i
0.0018 - 0.0123i 0.1711 + 0.0148i
0.0252 - 0.0173i 0.0428 + 0.0726i
0.0116 + 0.0120i -0.2852 + 0.0160i
-0.0390 + 0.0386i -0.2459 - 0.1224i
-0.0441 + 0.0009i 0.3736 - 0.1013i
0.0450 - 0.0697i 0.6682 + 0.1631i
0.1069 - 0.0427i -0.2866 + 0.2809i
-0.0190 + 0.1019i -1.3887 - 0.1317i
-0.2065 + 0.1383i -0.2970 - 0.5910i
-0.0874 - 0.1046i 2.3639 - 0.1084i
params = lsqcurvefit(@(params, r) [dr(r, params(1),params(2),params(3),params(4)),dtheta(r,params(1),params(2),params(3),params(4))], [lambda_init, kappa_init,theta_k_init,R_init], [dr_data(:, 1),dtheta_data(:,1)], [dr_data(:, 2),dtheta_data(:,2)], lb, ub, options);
Error using lsqncommon (line 65)
Lower and upper bounds not supported with complex-valued initial function or Jacobian evaluation.
Lower and upper bounds not supported with complex-valued initial function or Jacobian evaluation.
Error in lsqcurvefit (line 306)
lsqncommon(funfcn,xCurrent,lb,ub,options,defaultopt,optimgetFlag,caller,...
% Extract the optimized parameters
lambda_fit = params(1);
kappa_fit = params(2);
theta_k_fit = params(3);
R_fit = params(4);
% Use the initial guess for theta_k and R
%theta_k_fit = theta_k;
%R_fit = R;
% Plotting data
r_values = linspace(min(dr_data(:, 1)), max(dr_data(:, 1)), 1000).';
dr_fit = dr([r_values,r_values], params(1),params(2),params(3),params(4));
dtheta_fit = dtheta([r_values,r_values], params(1),params(2),params(3),params(4));
% Plotting
figure;
subplot(2, 1, 1);
plot(dr_data(:, 1), dr_data(:, 2), 'bo', 'DisplayName', 'Data');
hold on;
plot(r_values, dr_fit, 'r-', 'DisplayName', 'Fitted');
xlabel('r');
ylabel('dr');
title('Fitting dr');
legend;
subplot(2, 1, 2);
plot(dtheta_data(:, 1), dtheta_data(:, 2), 'bo', 'DisplayName', 'Data');
hold on;
plot(r_values, dtheta_fit, 'r-', 'DisplayName', 'Fitted');
xlabel('r');
ylabel('dtheta');
title('Fitting dtheta');
legend;
Hi @Sam,
Thanks for comment.
I could not able to fix it. There is nothing wrong in the function.
It's some technical error. you can check the error messege.
Torsten
am 2 Apr. 2024
⋮
-
Verknüpfen
Direkter Link zu diesem Kommentar
Bearbeitet: Torsten
am 2 Apr. 2024
The critical expression seems to be
(kappa^2 - omega_m(lambda, kappa, theta_k)^4) * ...
(kappa^2 - omega_p(lambda, kappa, theta_k)^4)
which has to be >= 0.
Arrange your initial parameters kappa_init,... such that this is the case.
Another possible source of complex numbers can be omega_m and omega_p where also sqrt expressions exist.
After adjusting the initial parameter values, use "nonlcon" to define the conditions as inequality constraints
to keep your function values real.
tuhin
am 2 Apr. 2024
⋮
-
Verknüpfen
Direkter Link zu diesem Kommentar
can you please modify the code and include this within the code?
Torsten
am 2 Apr. 2024
⋮
-
Verknüpfen
Direkter Link zu diesem Kommentar
No, sorry. This is too much work to do.
Siehe auch
Kategorien
Mehr zu Interpolation of 2-D Selections in 3-D Grids finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!Es ist ein Fehler aufgetreten
Da Änderungen an der Seite vorgenommen wurden, kann diese Aktion nicht abgeschlossen werden. Laden Sie die Seite neu, um sie im aktualisierten Zustand anzuzeigen.
Translated by 
Website auswählen
Wählen Sie eine Website aus, um übersetzte Inhalte (sofern verfügbar) sowie lokale Veranstaltungen und Angebote anzuzeigen. Auf der Grundlage Ihres Standorts empfehlen wir Ihnen die folgende Auswahl: .
Sie können auch eine Website aus der folgenden Liste auswählen:
So erhalten Sie die bestmögliche Leistung auf der Website
Wählen Sie für die bestmögliche Website-Leistung die Website für China (auf Chinesisch oder Englisch). Andere landesspezifische Websites von MathWorks sind für Besuche von Ihrem Standort aus nicht optimiert.
Amerika
- América Latina (Español)
- Canada (English)
- United States (English)
Europa
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
Asien-Pazifik
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)