Fitting not working or bad fit

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

Answer 1

[akshatsood]

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.