You must fit both data curves simultaneously, not in two separate calls to "lsqcurvefit".
You are now following this question
- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.
Simulation data Fittting problem
1 view (last 30 days)
Show older comments
% 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 Comments
tuhin
on 1 Apr 2024
Edited: Voss
on 1 Apr 2024
% Rearranged Data Matrix
dr_data = [2.1714, -0.0059213; 4.3429, 0.017543; 6.5143, 0.086; 8.6857, 0.15588; 10.8571, 0.20539; 13.0286, 0.19694; 15.2, 0.20442; 17.3714, 0.15659; 19.5429, 0.18918; 21.7143, 0.18262; 23.8857, 0.13818; 26.0571, 0.1539; 28.2286, 0.11195; 30.4, 0.12689; 32.5714, 0.068146; 34.7429, 0.028182; 36.9143, 0.013852; 39.0857, 0.039137; 41.2571, 0.00033664; 43.4286, -0.036782; 45.6, -0.043573; 47.7714, -0.060933; 49.9429, -0.030135; 52.1143, -0.043654; 54.2857, -0.039393; 56.4571, -0.030637; 58.6286, -0.03931; 60.8, -0.044883; 62.9714, -0.022349; 65.1429, -0.01046; 67.3143, 0.0014764; 69.4857, 0.012712; 71.6571, 0.0211; 73.8286, 0.02204];
dtheta_data = [2.1714, -0.011123; 4.3429, -0.31772; 6.5143, -0.3745; 8.6857, -0.40013; 10.8571, -0.50617; 13.0286, -0.49345; 15.2, -0.44292; 17.3714, -0.42858; 19.5429, -0.41354; 21.7143, -0.29636; 23.8857, -0.22671; 26.0571, -0.099143; 28.2286, 0.0087113; 30.4, -0.042737; 32.5714, 0.01474; 34.7429, 0.11353; 36.9143, 0.094084; 39.0857, 0.11759; 41.2571, 0.16252; 43.4286, 0.16718; 45.6, 0.22171; 47.7714, 0.25543; 49.9429, 0.25836; 52.1143, 0.23052; 54.2857, 0.12648; 56.4571, 0.15518; 58.6286, 0.20362; 60.8, 0.22967; 62.9714, 0.22242; 65.1429, 0.19043; 67.3143, 0.1749; 69.4857, 0.16304; 71.6571, 0.14256; 73.8286, 0.14299];
% 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
Answers (1)
Torsten
on 1 Apr 2024
Edited: Torsten
on 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 Comments
tuhin
on 1 Apr 2024
% Rearranged Data Matrix
dr_data = [2.1714, -0.0059213; 4.3429, 0.017543; 6.5143, 0.086; 8.6857, 0.15588; 10.8571, 0.20539; 13.0286, 0.19694; 15.2, 0.20442; 17.3714, 0.15659; 19.5429, 0.18918; 21.7143, 0.18262; 23.8857, 0.13818; 26.0571, 0.1539; 28.2286, 0.11195; 30.4, 0.12689; 32.5714, 0.068146; 34.7429, 0.028182; 36.9143, 0.013852; 39.0857, 0.039137; 41.2571, 0.00033664; 43.4286, -0.036782; 45.6, -0.043573; 47.7714, -0.060933; 49.9429, -0.030135; 52.1143, -0.043654; 54.2857, -0.039393; 56.4571, -0.030637; 58.6286, -0.03931; 60.8, -0.044883; 62.9714, -0.022349; 65.1429, -0.01046; 67.3143, 0.0014764; 69.4857, 0.012712; 71.6571, 0.0211; 73.8286, 0.02204];
dtheta_data = [2.1714, -0.011123; 4.3429, -0.31772; 6.5143, -0.3745; 8.6857, -0.40013; 10.8571, -0.50617; 13.0286, -0.49345; 15.2, -0.44292; 17.3714, -0.42858; 19.5429, -0.41354; 21.7143, -0.29636; 23.8857, -0.22671; 26.0571, -0.099143; 28.2286, 0.0087113; 30.4, -0.042737; 32.5714, 0.01474; 34.7429, 0.11353; 36.9143, 0.094084; 39.0857, 0.11759; 41.2571, 0.16252; 43.4286, 0.16718; 45.6, 0.22171; 47.7714, 0.25543; 49.9429, 0.25836; 52.1143, 0.23052; 54.2857, 0.12648; 56.4571, 0.15518; 58.6286, 0.20362; 60.8, 0.22967; 62.9714, 0.22242; 65.1429, 0.19043; 67.3143, 0.1749; 69.4857, 0.16304; 71.6571, 0.14256; 73.8286, 0.14299];
% 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
on 1 Apr 2024
Edited: Torsten
on 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
on 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
on 1 Apr 2024
Edited: tuhin
on 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
on 1 Apr 2024
Edited: Torsten
on 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
on 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
on 2 Apr 2024
Edited: Torsten
on 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
on 2 Apr 2024
Edited: tuhin
on 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
on 2 Apr 2024
Edited: Torsten
on 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
on 2 Apr 2024
Edited: Torsten
on 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
on 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
on 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
on 2 Apr 2024
Edited: Torsten
on 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.
See Also
Categories
Find more on Parallel Computing Fundamentals in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!An Error Occurred
Unable to complete the action because of changes made to the page. Reload the page to see its updated state.
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list
How to Get Best Site Performance
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.
Americas
- 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)
Asia Pacific
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)