solving system of 4 coupled odes using shooting method and boundary conditions given
4 views (last 30 days)
Show older comments
Hello! i am trying to solve a boundary value problem with four coupled first order odes, with four initial conditions at r=0 and four boundary conditions at r=10. my code is running but im not getting my desired output and it is not satisfying the conditions. Can someone help me finding where am i going wrong?
function proca_star_shooting_method
% Clear the workspace and command window
clear;
clc;
format long;
% Define constants and parameters
infinity = 10;
w = 0.817;
x_init = linspace(1e-5, infinity, 1000); % More efficient space vector
% Define initial conditions and boundary conditions
initial_conditions = [0.394, 0.394, 0, 0];
boundary_conditions = [1, 0, 0.745, 0];
% Shooting method with RK4 integration
options = optimset('TolX', 1e-6, 'Display', 'iter');
phi_c_shoot = fminbnd(@(phi_c) shooting_error(phi_c, x_init, initial_conditions, boundary_conditions, w), 0, 1, options);
% Solve the BVP using the optimal phi_c obtained from the shooting method
[r_data, y_data] = solve_bvp(x_init, initial_conditions, phi_c_shoot, w);
% Plot the results
plot_results(r_data, y_data, infinity);
end
function error = shooting_error(phi_c, x_init, initial_conditions, boundary_conditions, w)
% Solve BVP with given phi_c and calculate the error at the boundary
[~, y_data] = solve_bvp(x_init, initial_conditions, phi_c, w);
yb = y_data(:, end);
error = norm(yb - boundary_conditions');
end
function [r_data, y_data] = solve_bvp(x_init, initial_conditions, phi_c, w)
% Define the initial conditions for the solver
y0 = [initial_conditions(1), initial_conditions(2), initial_conditions(3), phi_c];
% Initialize the solution arrays
r_data = x_init;
y_data = zeros(4, length(x_init));
y_data(:, 1) = y0;
h = x_init(2) - x_init(1); % Step size
% RK4 integration loop
for i = 2:length(x_init)
y_data(:, i) = rk4_step(@(r, y) bsode(r, y, w), r_data(i-1), y_data(:, i-1), h);
end
end
function y_next = rk4_step(odefun, r, y, h)
% Perform one step of RK4 integration
k1 = h * odefun(r, y);
k2 = h * odefun(r + h/2, y + k1/2);
k3 = h * odefun(r + h/2, y + k2/2);
k4 = h * odefun(r + h, y + k3);
y_next = y + (k1 + 2*k2 + 2*k3 + k4) / 6;
end
function dfdr = bsode(r, y, w)
% Define the system of ODEs to be solved
N = 1 - 2 * y(3) / r;
dfdr = [
4 * pi * r * 0.745 * (y(4)^2 + y(2)^2 / (N^2 * y(1)^2));
w * y(4) - 0.745 * y(1)^2 * N * y(4) / w;
4 * pi * r^2 * (((0.745 * y(1)^2 * N^2 * y(4)^2)^2) / (2 * y(1)^2 * w^2) + 0.5 * 0.745 * (y(4)^2 * N + y(2)^2 / (N * y(1)^2)));
r^2 * y(2) * w / (y(1)^2 * N^2) - 2 * y(4)
];
end
function plot_results(r_data, y_data, infinity)
% Plot the results of the BVP solution
figure;
plot(r_data, y_data(1,:), r_data, y_data(2,:), r_data, y_data(3,:), r_data, y_data(4,:));
axis([0 infinity -0.5 1.2]);
title('Functions vs r');
xlabel('r');
ylabel('\sigma');
legend('y1', 'y2', 'y3', 'y4');
grid on;
end
This is my desired plot. Please help me out!![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1735462/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1735462/image.png)
1 Comment
J. Alex Lee
on 17 Jul 2024
just by counting conditions, if you only have 4 x 1st order ODEs, you can only have 4 BCs, you can't have 8?
Answers (0)
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!