Problem with the fitting

% 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 initial guess values for parameters
lambda_init = 2.9;
kappa_init = 0.066;
theta_k_init = 0.45;
R_init = 2400;
rout = 76.3;

% Define c
c = sqrt(4 - (rout/R_init)^2);

% Define initial conditions as parameters
a0_bound = 0.0;
b0_bound = 0.0;
c0_bound = 0.0;
d0_bound = 0.0;

% Objective function
objective_function = @(params) compute_error(params, dr_data, dtheta_data, R_init, rout, c);

% Optimization options
options = optimoptions(@fminunc,'Display','iter','Algorithm','quasi-newton','MaxFunctionEvaluations',10000,'MaxIterations',10000);

% Perform optimization
best_params = fminunc(objective_function, [lambda_init, kappa_init, theta_k_init, a0_bound, b0_bound, c0_bound, d0_bound], options);

% Display best parameters
disp("Optimal Parameters:");
disp("Lambda: " + best_params(1));
disp("Kappa: " + best_params(2));
disp("Theta_k: " + best_params(3));

% Plot the results with the best parameters
[~, dr_fit, dtheta_fit] = compute_error(best_params, dr_data, dtheta_data, R_init, rout, c);

figure;
subplot(2,1,1);
plot(dr_data(:,1), dr_data(:,2), 'ro', dr_data(:,1), dr_fit, 'b-');
xlabel('r');
ylabel('dr(r)');
title('Fit of dr(r) vs r');
legend('Data', 'Fit');

subplot(2,1,2);
plot(dtheta_data(:,1), dtheta_data(:,2), 'ro', dtheta_data(:,1), dtheta_fit, 'b-');
xlabel('r');
ylabel('dtheta(r)');
title('Fit of dtheta(r) vs r');
legend('Data', 'Fit');

% Error computation function
function [error, dr_fit, dtheta_fit] = compute_error(params, dr_data, dtheta_data, R, rout, c)
    lambda = params(1);
    kappa = params(2);
    theta_k = params(3);
    a0 = params(4);
    b0 = params(5);
    c0 = params(6);
    d0 = params(7);

    % Define the function for the differential equations
    f = @(r, y) [y(2); ((-1*(lambda+1)*(r*y(2)+y(1)))+(1/r)*((kappa*r^2*cos(theta_k)-(lambda+1))*y(1))-(kappa*r*y(3)*sin(theta_k))+(-16*lambda*r^2)/(c^4*R^2)); y(4); ((-1*(r*y(4)+y(3)))+(1/r)*((kappa*r^2*cos(theta_k)-1)*y(3))+(kappa*r*y(1)*sin(theta_k)))];

    % Solve the differential equations using ode45
    %r_span = [min(dr_data(:,1)), max(dr_data(:,1))]; % Span of r values from the data
    r_span = dr_data(:,1);
    [~, sol] = ode45(f, r_span, [a0, b0, c0, d0]);

    % Extract solutions
    dr_fit = sol(:,1);
    dtheta_fit = sol(:,3);
    
    % Compute error (sum of squared differences)
    error_dr = sum((dr_data(:,2) - dr_fit).^2);
    error_dtheta = sum((dtheta_data(:,2) - dtheta_fit).^2);
    
    % Total error
    error = error_dr + error_dtheta;
end

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 also. However getting errors. Please help me to solve those errors and fitting those data.

The mathematical eqns and details are below:

$\kappa>0$ and $0 \leq \theta_{k}<\frac{\pi}{2}$ describe the screening. The boundary conditions are $\vec{d}(0)=\vec{d}\left(r_{\text {out }}\right)=(0,0)$. 

For 
$$
c=\sqrt{4-\left(\frac{r_{\text {out }}}{R}\right)^{2}}
$$

Furthermore, 

\begin{equation}
(\lambda+1)\left(r d_{r}^{\prime \prime}(r)+d_{r}^{\prime}(r)\right) &+ \frac{1}{r}\left(\kappa r^{2} \cos \theta_{k}-(\lambda+1)\right) d_{r}(r)\\
& - \kappa r d_{\theta}(r) \sin \theta_{k} = -\frac{16 \lambda r^{2}}{c^{4} R^{2}}
\end{equation}

\begin{equation}
r d_{\theta}^{\prime \prime}(r) + d_{\theta}^{\prime}(r) &+ \frac{1}{r}\left(\kappa r^{2} \cos \theta_{k}-1\right) d_{\theta}(r)\\
& + \kappa r d_{r}(r) \sin \theta_{k} = 0
\end{equation}

[Adam Danz]

I ran your code to produce the error message. Is this the error you receive? Arrays have incompatible sizes for this operation.

Answer 1

[Torsten]

% 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 initial guess values for parameters
lambda_init = 2.9;
kappa_init = 0.066;
theta_k_init = 0.45;
R_init = 2400;
rout = 76.3;

% Define c
c = sqrt(4 - (rout/R_init)^2);

% Define initial conditions as parameters
a0_bound = 0.0;
b0_bound = 0.0;
c0_bound = 0.0;
d0_bound = 0.0;

% Objective function
objective_function = @(params) compute_error(params, dr_data, dtheta_data, R_init, rout, c);

% Optimization options
options = optimoptions(@fminunc,'Display','iter','Algorithm','quasi-newton','MaxFunctionEvaluations',10000,'MaxIterations',10000);

% Perform optimization
best_params = fminunc(objective_function, [lambda_init, kappa_init, theta_k_init, a0_bound, b0_bound, c0_bound, d0_bound], options);

                                                        First-order 
 Iteration  Func-count       f(x)        Step-size       optimality
     0           8         0.661751                          1.47
     1          32         0.603093      0.0268002           1.19  
     2          56         0.601328       0.111336            0.5  
     3          64         0.596533              1          0.394  
     4          80         0.588062            0.5           0.53  
     5          88         0.584744              1          0.209  
     6          96         0.582801              1           1.02  
     7         104         0.580027              1          0.914  
     8         112         0.574907              1          0.536  
     9         120         0.572701              1          0.739  
    10         128          0.57137              1          0.254  
    11         136         0.568637              1          0.481  
    12         168         0.510233         4.9344           2.11  
    13         192         0.505129      0.0334871           2.47  
    14         200         0.477892              1           2.22  
    15         208         0.451203              1            2.5  
    16         224         0.405833       0.493954           2.11  
    17         240         0.393151            0.1           1.03  
    18         248         0.374132              1           0.73  
    19         264         0.353052            0.5          0.888  
                                                        First-order 
 Iteration  Func-count       f(x)        Step-size       optimality
    20         280         0.349712       0.388692          0.913  
    21         288          0.34634              1          0.244  
    22         296         0.344831              1          0.213  
    23         304         0.344082              1          0.591  
    24         312         0.342182              1          0.873  
    25         320         0.338108              1           1.36  
    26         328         0.331145              1            1.7  
    27         336         0.323695              1           1.59  
    28         344         0.311484              1          0.495  
    29         352         0.295257              1          0.488  
    30         360         0.266355              1           1.14  
    31         368         0.252549              1           1.69  
    32         376          0.24712              1           0.66  
    33         384         0.244712              1          0.254  
    34         392         0.243488              1          0.276  
    35         400         0.243275              1          0.168  
    36         408         0.243127              1          0.226  
    37         416         0.242744              1          0.264  
    38         424         0.242058              1          0.279  
    39         432         0.241334              1          0.204  
                                                        First-order 
 Iteration  Func-count       f(x)        Step-size       optimality
    40         440         0.240979              1         0.0733  
    41         448         0.240919              1          0.013  
    42         456          0.24091              1         0.0131  
    43         464         0.240896              1         0.0227  
    44         472          0.24086              1         0.0407  
    45         480         0.240769              1         0.0688  
    46         488          0.24053              1          0.113  
    47         496         0.239929              1          0.178  
    48         504         0.238495              1          0.259  
    49         512         0.235441              1          0.343  
    50         520         0.229949              1          0.444  
    51         528         0.222181              1          0.624  
    52         536         0.215689              1          0.634  
    53         544         0.213293              1          0.303  
    54         552         0.212665              1          0.247  
    55         560         0.212384              1          0.024  
    56         568         0.212322              1         0.0291  
    57         576         0.212312              1         0.0212  
    58         584         0.212299              1          0.014  
    59         592         0.212298              1        0.00637  
                                                        First-order 
 Iteration  Func-count       f(x)        Step-size       optimality
    60         600         0.212298              1        0.00317  
    61         608         0.212297              1        0.00697  
    62         616         0.212294              1         0.0184  
    63         624         0.212287              1         0.0384  
    64         632         0.212266              1         0.0717  
    65         640         0.212222              1          0.114  
    66         648         0.212114              1          0.175  
    67         656         0.211833              1           0.28  
    68         664         0.211304              1          0.382  
    69         672         0.210304              1          0.441  
    70         680         0.209126              1          0.318  
    71         688         0.208632              1          0.311  
    72         704         0.208485       0.325369         0.0849  
    73         712           0.2084              1         0.0658  
    74         720         0.208228              1         0.0381  
    75         728         0.208173              1         0.0272  
    76         736         0.208134              1          0.022  
    77         744         0.208126              1         0.0033  
    78         752         0.208125              1         0.0118  
    79         760         0.208125              1       0.000156  
                                                        First-order 
 Iteration  Func-count       f(x)        Step-size       optimality
    80         768         0.208125              1        0.00018  
    81         784         0.208125       0.298081       0.000181  
    82         792         0.208125              1       8.68e-05  
    83         800         0.208125              1       1.43e-05  
    84         808         0.208125              1       5.72e-06  
    85         824         0.208125       0.396005       2.86e-06  
    86         832         0.208125              1       6.89e-08  

Local minimum found.

Optimization completed because the size of the gradient is less than
the value of the optimality tolerance.

% Display best parameters
disp("Optimal Parameters:");

Optimal Parameters:

disp("Lambda: " + best_params(1));

Lambda: 11.6368

disp("Kappa: " + best_params(2));

Kappa: 0.26924

disp("Theta_k: " + best_params(3));

Theta_k: 1.2933

% Plot the results with the best parameters
[~, dr_fit, dtheta_fit] = compute_error(best_params, dr_data, dtheta_data, R_init, rout, c);

figure;
subplot(2,1,1);
plot(dr_data(:,1), dr_data(:,2), 'ro', dr_data(:,1), dr_fit, 'b-');
xlabel('r');
ylabel('dr(r)');
title('Fit of dr(r) vs r');
legend('Data', 'Fit');

subplot(2,1,2);
plot(dtheta_data(:,1), dtheta_data(:,2), 'ro', dtheta_data(:,1), dtheta_fit, 'b-');
xlabel('r');
ylabel('dtheta(r)');
title('Fit of dtheta(r) vs r');
legend('Data', 'Fit');

% Error computation function
function [error, dr_fit, dtheta_fit] = compute_error(params, dr_data, dtheta_data, R, rout, c)
    lambda = params(1);
    kappa = params(2);
    theta_k = params(3);
    a0 = params(4);
    b0 = params(5);
    c0 = params(6);
    d0 = params(7);

    % Define the function for the differential equations
    f = @(r, y) [y(2); ((-1*(lambda+1)*(r*y(2)+y(1)))+(1/r)*((kappa*r^2*cos(theta_k)-(lambda+1))*y(1))-(kappa*r*y(3)*sin(theta_k))+(-16*lambda*r^2)/(c^4*R^2)); y(4); ((-1*(r*y(4)+y(3)))+(1/r)*((kappa*r^2*cos(theta_k)-1)*y(3))+(kappa*r*y(1)*sin(theta_k)))];

    % Solve the differential equations using ode45
    %r_span = [min(dr_data(:,1)), max(dr_data(:,1))]; % Span of r values from the data
    r_span = dr_data(:,1);
    [~, sol] = ode45(f, r_span, [a0, b0, c0, d0]);

    % Extract solutions
    dr_fit = sol(:,1);
    dtheta_fit = sol(:,3);
    
    % Compute error (sum of squared differences)
    error_dr = sum((dr_data(:,2) - dr_fit).^2);
    error_dtheta = sum((dtheta_data(:,2) - dtheta_fit).^2);
    
    % Total error
    error = error_dr + error_dtheta;
end