data not generating and error in plotting

1 次查看（过去 30 天）
tuhin
tuhin 2024-5-30
编辑： Torsten 2024-5-30
Ran in:
% 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.

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

采纳的回答

Torsten
Torsten 2024-5-30
移动：Torsten 2024-5-30
Ran in:
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
  2 个评论
tuhin
tuhin 2024-5-30
编辑：Torsten 2024-5-30
Ran in:
% 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
Torsten 2024-5-30
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.

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Line Plots 的更多信息

标签

产品


版本

R2023b

Community Treasure Hunt

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

Start Hunting!

Translated by