You must fit both data curves simultaneously, not in two separate calls to "lsqcurvefit".
Al momento, stai seguendo questa domanda
- Vedrai gli aggiornamenti nel tuofeed del contenuto seguito.
- Potresti ricevere delle e-mail a seconda delle tuepreferenze per le comunicazioni.
Simulation data Fittting problem
1 visualizzazione (ultimi 30 giorni)
Mostra commenti meno recenti
% 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 Commenti
tuhin
il 1 Apr 2024
Modificato: Voss
il 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
Risposte (1)
Torsten
il 1 Apr 2024
Modificato: Torsten
il 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 Commenti
Torsten
il 1 Apr 2024
Modificato: Torsten
il 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
il 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];
% 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
il 1 Apr 2024
Modificato: Torsten
il 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
il 1 Apr 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
il 1 Apr 2024
Modificato: tuhin
il 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
il 1 Apr 2024
Modificato: Torsten
il 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
il 1 Apr 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
il 2 Apr 2024
Modificato: Torsten
il 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
il 2 Apr 2024
Modificato: tuhin
il 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
il 2 Apr 2024
Modificato: Torsten
il 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
il 2 Apr 2024
Modificato: Torsten
il 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.
Sam Chak
il 2 Apr 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
il 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_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
il 2 Apr 2024
Modificato: Torsten
il 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.
Vedere anche
Categorie
Scopri di più su Parallel Computing Fundamentals in Help Center e File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!Si è verificato un errore
Impossibile completare l'azione a causa delle modifiche apportate alla pagina. Ricarica la pagina per vedere lo stato aggiornato.
Translated by 
Seleziona un sito web
Seleziona un sito web per visualizzare contenuto tradotto dove disponibile e vedere eventi e offerte locali. In base alla tua area geografica, ti consigliamo di selezionare: .
Puoi anche selezionare un sito web dal seguente elenco:
Come ottenere le migliori prestazioni del sito
Per ottenere le migliori prestazioni del sito, seleziona il sito cinese (in cinese o in inglese). I siti MathWorks per gli altri paesi non sono ottimizzati per essere visitati dalla tua area geografica.
Americhe
- 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)
Asia-Pacifico
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)