You must fit both data curves simultaneously, not in two separate calls to "lsqcurvefit".
Vous suivez désormais cette question
- Les mises à jour seront visibles dans votre flux de contenu suivi.
- Selon vos préférences en matière de communication il est possible que vous receviez des e-mails.
Simulation data Fittting problem
1 vue (au cours des 30 derniers jours)
Afficher commentaires plus anciens
% 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 commentaires
tuhin
le 1 Avr 2024
Modifié(e) : Voss
le 1 Avr 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
Réponses (1)
Torsten
le 1 Avr 2024
Modifié(e) : Torsten
le 1 Avr 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 commentaires
Torsten
le 1 Avr 2024
Modifié(e) : Torsten
le 1 Avr 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
le 1 Avr 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];
% 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
le 1 Avr 2024
Modifié(e) : Torsten
le 1 Avr 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
le 1 Avr 2024
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
le 1 Avr 2024
Modifié(e) : tuhin
le 1 Avr 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
le 1 Avr 2024
Modifié(e) : Torsten
le 1 Avr 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
le 1 Avr 2024
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
le 2 Avr 2024
Modifié(e) : Torsten
le 2 Avr 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
le 2 Avr 2024
Modifié(e) : tuhin
le 2 Avr 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
le 2 Avr 2024
Modifié(e) : Torsten
le 2 Avr 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
le 2 Avr 2024
Modifié(e) : Torsten
le 2 Avr 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.
Sam Chak
le 2 Avr 2024
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
le 2 Avr 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_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
le 2 Avr 2024
Modifié(e) : Torsten
le 2 Avr 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.
Voir également
Catégories
En savoir plus sur Parallel Computing Fundamentals dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!Une erreur s'est produite
Impossible de terminer l’action en raison de modifications de la page. Rechargez la page pour voir sa mise à jour.
Translated by 
Sélectionner un site web
Choisissez un site web pour accéder au contenu traduit dans votre langue (lorsqu'il est disponible) et voir les événements et les offres locales. D’après votre position, nous vous recommandons de sélectionner la région suivante : .
Vous pouvez également sélectionner un site web dans la liste suivante :
Comment optimiser les performances du site
Pour optimiser les performances du site, sélectionnez la région Chine (en chinois ou en anglais). Les sites de MathWorks pour les autres pays ne sont pas optimisés pour les visites provenant de votre région.
Amériques
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- 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)
Asie-Pacifique
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)