Fitting not working or bad fit

2 次查看（过去 30 天）
tuhin
tuhin 2024-5-30
编辑： akshatsood 2024-6-22
% Define the data
dr_data = [2.5453, 0.042123; 5.0907, 0.075326; 7.636, 0.059506; 10.1813, 0.071553; 12.7267, 0.071365; 15.272, 0.067195; 17.8173, 0.046372; 20.3627, 0.043397; 22.908, 0.017179; 25.4533, -0.0063329; 27.9987, -0.030789; 30.544, -0.047569; 33.0893, -0.089512; 35.6347, -0.080675; 38.18, -0.089138; 40.7253, -0.1102; 43.2707, -0.12061; 45.816, -0.11857; 48.3613, -0.11955; 50.9067, -0.10803; 53.452, -0.10462; 55.9973, -0.099548; 58.5427, -0.097164; 61.088, -0.09994; 63.6333, -0.077017; 66.1787, -0.062839; 68.724, -0.048422; 71.2693, -0.03686; 73.8147, -0.01469; 76.3, 0];
dtheta_data = [2.5453, -0.099251; 5.0907, -0.16064; 7.636, -0.21858; 10.1813, -0.18965; 12.7267, -0.16996; 15.272, -0.18172; 17.8173, -0.15029; 20.3627, -0.12541; 22.908, -0.082786; 25.4533, -0.0071716; 27.9987, 0.03695; 30.544, 0.089002; 33.0893, 0.12873; 35.6347, 0.13092; 38.18, 0.13908; 40.7253, 0.17211; 43.2707, 0.16686; 45.816, 0.15826; 48.3613, 0.14872; 50.9067, 0.15295; 53.452, 0.12677; 55.9973, 0.10964; 58.5427, 0.10223; 61.088, 0.10951; 63.6333, 0.088493; 66.1787, 0.068903; 68.724, 0.054396; 71.2693, 0.035731; 73.8147, 0.030172; 76.3, 0];
% Define the function for the differential equations
function dydr = odefun(r, y, lambda_init, kappa_init, theta_k_init, c, R)
dydr = zeros(4,1);
dydr(1) = y(2);
dydr(2) = -((lambda_init+1)*y(2)+1/r*(kappa_init*r^2*cos(theta_k_init)-(lambda_init+1))*y(1)-kappa_init*r*y(3)*sin(theta_k_init)+16*lambda_init*r^2/(c^4*R^2))/(r*(lambda_init+1));
dydr(3) = y(4);
dydr(4) = -(y(4)+1/r*(kappa_init*r^2*cos(theta_k_init)-1)*y(3)+kappa_init*r*y(1)*sin(theta_k_init))/r;
end
% Boundary conditions
function res = bcfun(ya, yb, ya1, ya3, yb1, yb3)
res = [ya(1)-ya1; ya(3)-ya3; yb(1)-yb1; yb(3)-yb3];
end
% Define the function to compute residuals
function residuals = compute_residuals(params, dr_data, dtheta_data, rout, R_init)
lambda_init = params(1);
kappa_init = params(2);
theta_k_init = params(3);
ya1 = params(4); % Boundary condition parameter ya(1)
ya3 = params(5); % Boundary condition parameter ya(3)
yb1 = params(6); % Boundary condition parameter yb(1)
yb3 = params(7); % Boundary condition parameter
c = sqrt(4 - (rout/R_init)^2);
R_init = 2000; % Initial value of R
rout = 76.3; % Max value of r
% Adjust the number of points for interpolation
num_points = size(dr_data, 1);
solinit = bvpinit(linspace(0.0001, rout, num_points), [ya1, 0, ya3, 0]);
sol = bvp4c(@(r, y) odefun(r, y, lambda_init, kappa_init, theta_k_init, c, R_init), @(ya, yb) bcfun(ya, yb, ya1, ya3, yb1, yb3), solinit);
r = linspace(0.0001, rout, num_points);
y = deval(sol, r);
dr_sol = y(1,:);
dtheta_sol = y(3,:);
% Ensure that arrays have compatible sizes
dr_residuals = dr_data(:, 2) - dr_sol;
dtheta_residuals = dtheta_data(:, 2) - dtheta_sol;
residuals = [dr_residuals; dtheta_residuals];
end
% Initial guess for the parameters
params_init = [2.1, 0.03, 0.004, 0, 0, 0, 0]; % Initial guess for parameters including boundary conditions
% Bounds for parameters
lb = [1, 0.001, 0, -Inf, -Inf, -Inf, -Inf]; % Lower bounds for lambda, kappa, theta_k, ya1, ya3, yb1, yb3
ub = [3, 0.5, pi/2, Inf, Inf, Inf, Inf]; % Upper bounds for lambda, kappa, theta_k, ya1, ya3, yb1, yb3
% Perform optimization
params_opt = lsqnonlin(@(params) compute_residuals(params, dr_data, dtheta_data, rout, R_init), params_init, lb, ub);
% Extract optimized parameters
lambda_opt = params_opt(1);
kappa_opt = params_opt(2);
theta_k_opt = params_opt(3);
ya1_opt = params_opt(4);
ya3_opt = params_opt(5);
yb1_opt = params_opt(6);
yb3_opt = params_opt(7);
% Display optimized parameters
disp(['Optimized lambda: ', num2str(lambda_opt)]);
disp(['Optimized kappa: ', num2str(kappa_opt)]);
disp(['Optimized theta_k: ', num2str(theta_k_opt)]);
disp(['Optimized ya1: ', num2str(ya1_opt)]);
disp(['Optimized ya3: ', num2str(ya3_opt)]);
disp(['Optimized yb1: ', num2str(yb1_opt)]);
disp(['Optimized yb3: ', num2str(yb3_opt)]);
% Plot the solutions using optimized parameters
lambda_init = lambda_opt;
kappa_init = kappa_opt;
theta_k_init = theta_k_opt;
ya1 = ya1_opt;
ya3 = ya3_opt;
yb1 = yb1_opt;
yb3 = yb3_opt;
c = sqrt(4 - (rout/R_init)^2);
num_points = size(dr_data, 1);
solinit = bvpinit(linspace(0.0001, rout, num_points), [ya1, yb1, ya3, yb3]);
sol = bvp4c(@(r, y) odefun(r, y, lambda_init, kappa_init, theta_k_init, c, R_init), @(ya, yb) bcfun(ya, yb, ya1, ya3, yb1, yb3), solinit);
r = linspace(0.0001, rout, num_points);
y = deval(sol, r);
dr_sol = y(1,:);
dtheta_sol = y(3,:);
% Plot the solutions
figure;
subplot(2,1,1);
plot(r, dr_sol, 'b-', dr_data(:,1), dr_data(:,2), 'ro');
xlabel('r');
ylabel('dr(r)');
title('Solution of dr(r) vs r');
legend('Fitted Solution', 'Data');
subplot(2,1,2);
plot(r, dtheta_sol, 'b-', dtheta_data(:,1), dtheta_data(:,2), 'ro');
xlabel('r');
ylabel('dtheta(r)');
title('Solution of dtheta(r) vs r');
legend('Fitted Solution', 'Data');
Now, I want to fit the simulated d_r(r) vs r and d_theta(r) vs r values with the above mentioned coupled differential eqns by usuing the fitting parameters lamda, kappa, theta_k. For the sake of good fitting one can use boundary conditions as parameter. However getting errors. Please help me to solve those errors and fitting those data.

请先登录，再回答此问题。

回答（1 个）

akshatsood
akshatsood 2024-6-22
编辑：akshatsood 2024-6-22
Dear @tuhin,
I understand that you are encountering errors while attempting to solve nonlinear least squares optimization problem.
Upon reviewing the code, I can correctly point out two major problems causing the errors:
  1. Scope of Variables "rout" and "R_init": The variables "rout" and "R_init" are defined inside the function "compute_residuals", making them inaccessible to other functions which require them as input.
  2. Function Definitions Placement: Function definitions should appear at the end of the file. The code that is not inside a function should be placed above the first function definition in your script.
Adhering to these two changes should be sufficient to execute the code without errors. However, I would recommend the following documentation for more information on solving nonlinear least squares problems.
https://www.mathworks.com/help/releases/R2023b/optim/nonlinear-least-squares-curve-fitting.html
Further, I have attached the updated script with this response. Kindly refer to it for the changes I have made.
I hope this helps.

类别

在 Help Center 和 File Exchange 中查找有关 Get Started with Optimization Toolbox 的更多信息

标签

产品


版本

R2023b

Community Treasure Hunt

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

Start Hunting!

Translated by