Help needed with fitting and complex number

% 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
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.

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,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;

Getting an error. Please help me to solve and fit it

Answer 1

[Nipun]

Hi Tuhin,

I understand that you are trying to fit the given data using the customized optimization functions but are facing an error with the "lsqcurvefit" function.

Based on the given error and the documentation on "lsqcurvefit" function, the function does not support complex valued functions with upper or lower bounds.

I suggest you either remove the bounds and try running the code to see if the output matches the expected output. Here is the MathWorks documentation on fitting complex data in MATLAB:

https://in.mathworks.com/help/optim/ug/fit-model-to-complex-data.html

Hope this helps.

Regards,

Nipun