Problem with the fitting
1 visualizzazione (ultimi 30 giorni)
Mostra commenti meno recenti
% 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:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1706506/image.jpeg)
$\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}
1 Commento
Adam Danz
il 30 Mag 2024
I ran your code to produce the error message. Is this the error you receive? Arrays have incompatible sizes for this operation.
Risposte (1)
Torsten
il 30 Mag 2024
% 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
0 Commenti
Vedere anche
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!