data not generating and error in plotting

% Define initial parameters
lambda_init = 2;
kappa_init = 1;
theta_k_init = pi/10;
R_init = 7;
rout = 3;

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

% Define the function for the differential equations
f = @(r, y) [y(2); ((-1*(lambda_init+1)*(r*y(2)+y(1)))+(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_init^2)); y(4); ((-1*(r*y(4)+y(3)))+(1/r)*((kappa_init*r^2*cos(theta_k_init)-1)*y(3))+(kappa_init*r*y(1)*sin(theta_k_init)))];

% Solve the differential equations using ode45
r_span = linspace(0, rout, 100); % Define the range of r values
[~, sol] = ode45(f, r_span, [0, 0, 0, 0]);

% Extract solutions
dr_sol = sol(:,1);
dtheta_sol = sol(:,3);

% Plot the solutions
figure;
subplot(2,1,1);
plot(r_span, dr_sol, 'b-');
xlabel('r');
ylabel('dr(r)');
title('Solution of dr(r) vs r');

subplot(2,1,2);
plot(r_span, dtheta_sol, 'b-');
xlabel('r');
ylabel('dtheta(r)');
title('Solution of dtheta(r) vs r');

I am trying to solve the coupled differential eqns for d_r(r) vs r and d_theta(r) vs r for these parameter values and boundary conditions so that it is zero at both end d_r(0) = d_r(rout) = 0 and the same for d_theta. However, Not able to see the plot and data. please suggest me errors.

Answer 1

[Torsten]

And what is the "correct" result ?

According to your mathematical description, I get this:

% Define initial parameters
lambda_init = 1.2;
kappa_init = 1;
theta_k_init = pi/10;
R_init = 7;
rout = 3;

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

% Initial guess for the solution
solinit = bvpinit(linspace(0.0001, rout, 100), [0, 0, 0, 0]);

% Solve the BVP
sol = bvp4c(@(r, y) odefun(r, y, lambda_init, kappa_init, theta_k_init, c, R_init), @bcfun, solinit);

% Extract solutions
r = linspace(0.0001, rout, 100);
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-');
xlabel('r');
ylabel('dr(r)');
title('Solution of dr(r) vs r');

subplot(2,1,2);
plot(r, dtheta_sol, 'b-');
xlabel('r');
ylabel('dtheta(r)');
title('Solution of dtheta(r) vs r');

% Define the function for the differential equations
function dydr = odefun(r, y, lambda_init, kappa_init, theta_k_init, c, R_init)
    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_init^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)
    res = [ya(1); ya(3); yb(1); yb(3)];
end

[tuhin]

% 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];

% Initial guess for the parameters
params_init = [1.2, 1, pi/10, 7];

% Bounds for parameters
lb = [1, 0.1, 0, 1]; % Lower bounds for lambda, kappa, theta_k, R
ub = [3, 10, pi/2, 10]; % Upper bounds for lambda, kappa, theta_k, R

% Define R_init and rout
R_init = 2000; % Initial value of R
rout = 76.3; % Max value of r

% Perform optimization
params_opt = lsqnonlin(@(params) compute_residuals(params, dr_data, dtheta_data, rout), params_init, lb, ub);

Local minimum possible.

lsqnonlin stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

% Extract optimized parameters
lambda_opt = params_opt(1);
kappa_opt = params_opt(2);
theta_k_opt = params_opt(3);
R_opt = params_opt(4);

% Display optimized parameters
disp(['Optimized lambda: ', num2str(lambda_opt)]);

Optimized lambda: 1.5586

disp(['Optimized kappa: ', num2str(kappa_opt)]);

Optimized kappa: 1.2095

disp(['Optimized theta_k: ', num2str(theta_k_opt)]);

Optimized theta_k: 0.38614

disp(['Optimized R: ', num2str(R_opt)]);

Optimized R: 9.0517

% Plot the solutions
lambda_init = lambda_opt;
kappa_init = kappa_opt;
theta_k_init = theta_k_opt;
R_init = R_opt;
c = sqrt(4 - (rout/R_init)^2);
solinit = bvpinit(linspace(0.0001, rout, 100), [0, 0, 0,0]);
sol = bvp4c(@(r, y) odefun(r, y, lambda_init, kappa_init, theta_k_init, c, R_init), @bcfun, solinit);
r = linspace(0.0001, rout, 100);
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');

% Define the function for the differential equations
function dydr = odefun(r, y, lambda_init, kappa_init, theta_k_init, c, R_init)
    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_init^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)
    res = [ya(1); ya(3); yb(1); yb(3)];
end

% Define the function to compute residuals
function residuals = compute_residuals(params, dr_data, dtheta_data, rout)
    lambda_init = params(1);
    kappa_init = params(2);
    theta_k_init = params(3);
    R_init = params(4);
    %rout = 76.3;
    c = sqrt(4 - (rout/R_init)^2);

    solinit = bvpinit(linspace(0.0001, rout, 100), [0, 0, 0, 0]);
    sol = bvp4c(@(r, y) odefun(r, y, lambda_init, kappa_init, theta_k_init, c, R_init), @bcfun, solinit);
    r = linspace(0.0001, rout, 100);
    y = deval(sol, dr_data(:,1));
    dr_sol = y(1,:);
    dtheta_sol = y(3,:);

    dr_residuals = dr_data(:, 2) - dr_sol.';
    dtheta_residuals = dtheta_data(:, 2) - dtheta_sol.';

    residuals = [dr_residuals; dtheta_residuals];
end

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.

[Torsten]

I changed

y = deval(sol, r);

to

y = deval(sol, dr_data(:,1));

to make your code work.

Fitting is bad - you'll have to work on it.