Help needed with fitting and complex number

3 次查看(过去 30 天)
% 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

回答(1 个)

Nipun
Nipun 2024-4-10
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:
Hope this helps.
Regards,
Nipun

类别

Help CenterFile Exchange 中查找有关 Computational Fluid Dynamics (CFD) 的更多信息

产品


版本

R2023b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by