getting singular jacobian error while using bvp4c solver
Show older comments
when i am using a range for phi_c for 0.01 to 0.1 , i am getting a plot. but when im increasing range from 0.01 to 0.8, i am getting error. please help me solve this.
% Define global variable for initial scalar field values
phi_c_values = logspace(log10(0.01), log10(0.1), 100);
% Physical constants
M_pl_squared = 1.0;
m = 1.0; % Scalar mass
% Function to define the differential equations
function dydr = bsode(r, y, p)
G = 1; % Gravitational constant
omega = p(1);
m=1;
a = y(1);
alpha = y(2);
psi0 = y(3);
phi = y(4);
dydr = zeros(4, 1);
dydr(1) = (0.5) * ((a/r) * ((1 - a^2) + 4 * pi * G * r * a * ...
(psi0^2 * a^2 * (m^2 + omega^2 / alpha^2) + phi^2)));
dydr(2) = (alpha/2) * (((a^2 - 1)/r) + 4 * pi * r * ...
(psi0^2 * a^2 * (omega^2 / alpha^2 - m^2) + phi^2));
dydr(3) = phi;
dydr(4) = -(1 + a^2 - 4 * pi * r^2 * psi0^2 * a^2 * m^2) * (phi/r) - ...
(omega^2 / alpha^2 - m^2) * psi0 * a^2;
end
% Function to define boundary conditions
function res = bsbc(ya, yb, p)
global phi_c;
res = [ya(1) - 1; % a(0) = 1
ya(2) - 1; % alpha(0) = 1
ya(3) - phi_c; % psi0(0) = phi_c
yb(3); % psi0(infinity) = 0
ya(4)]; % phi(0) = 0
end
% Function for calculating ADM mass
function M = ADM_mass(r, a)
M = (1 - a.^-2) .* (r / 2);
end
% Main loop over phi_c values
omega = 0.864;
infinity = 15;
x_init = linspace(1e-5, infinity, 1000);
% Lists to store maximum masses and radii
max_masses = [];
max_radii = [];
phi_values = [];
mass_values = [];
for phi_c = phi_c_values
% Initial guess structure
solinit = bvpinit(x_init, [1, 1, phi_c, 0], omega);
% Solve the BVP
sol = bvp4c(@bsode, @bsbc, solinit);
f = sol.y;
i = find(f(3, :) < 1e-5, 1);
r = sol.x(1:i);
a = f(1, 1:i);
% Scale alpha and omega
alpha_last = 1 / f(2, i);
scaled_alpha = alpha_last * f(2, 1:i);
scaled_omega = alpha_last * sol.parameters(1);
% Calculate ADM mass
M = ADM_mass(r, a);
M_max = 0.633 * M_pl_squared / m;
M_scaled = M / M_max;
% Find maximum mass and corresponding radius
max_mass = max(M_scaled);
target_mass = 0.99 * max_mass;
index = find(abs(M_scaled - target_mass) == min(abs(M_scaled - target_mass)), 1);
max_radius = r(index);
max_M = M_scaled(index);
% Store results
max_masses(end+1) = max_M;
max_radii(end+1) = max_radius;
phi_values(end+1) = phi_c;
mass_values(end+1) = max_mass;
fprintf('phi_c = %.5f, max_mass = %.5f, max_radius = %.5f\n', phi_c, max_M, max_radius);
end
% Plotting mass-radius curve
figure;
plot(max_radii, max_masses);
title('Scaled ADM Mass (M) vs Radial Coordinate (r)');
xlabel('r (Radius)');
ylabel('M (Scaled Mass)');
grid on;
% Plotting m vs phi
figure;
plot(phi_values, mass_values, 'orange');
title('Scaled ADM Mass (M) vs \phi_c');
xlabel('\phi_c');
ylabel('M (Scaled Mass)');
grid on;
Answers (1)
I can't interprete your solution curves, but it seems to be a problem for the solver that the Scaled Mass goes to 0 with increasing phi_c.
% Define global variable for initial scalar field values
phi_c_values = linspace(0.01,0.15,200);
% Physical constants
M_pl_squared = 1.0;
m = 1.0; % Scalar mass
% Main loop over phi_c values
omega = 0.864;
infinity = 15;
x_init = linspace(1e-5, infinity, 1000);
% Lists to store maximum masses and radii
max_masses = [];
max_radii = [];
phi_values = [];
mass_values = [];
solinit = bvpinit(x_init, [1, 1, phi_c_values(1), 0], omega);
options = bvpset('RelTol',1e-8,'AbsTol',1e-8);
for phi_c = phi_c_values
% Solve the BVP
sol = bvp4c(@bsode, @(x,y,p)bsbc(x,y,p,phi_c), solinit, options);
f = sol.y;
i = find(f(3, :) < 1e-5, 1);
r = sol.x(1:i);
a = f(1, 1:i);
% Scale alpha and omega
alpha_last = 1 / f(2, i);
scaled_alpha = alpha_last * f(2, 1:i);
scaled_omega = alpha_last * sol.parameters(1);
% Calculate ADM mass
M = ADM_mass(r, a);
M_max = 0.633 * M_pl_squared / m;
M_scaled = M / M_max;
% Find maximum mass and corresponding radius
max_mass = max(M_scaled);
target_mass = 0.99 * max_mass;
index = find(abs(M_scaled - target_mass) == min(abs(M_scaled - target_mass)), 1);
max_radius = r(index);
max_M = M_scaled(index);
% Store results
max_masses(end+1) = max_M;
max_radii(end+1) = max_radius;
phi_values(end+1) = phi_c;
mass_values(end+1) = max_mass;
fprintf('phi_c = %.5f, max_mass = %.5f, max_radius = %.5f\n', phi_c, max_M, max_radius);
solinit = bvpinit(sol.x,@(x)interp1(sol.x,sol.y.',x),sol.parameters(1));
end
% Plotting mass-radius curve
figure;
plot(max_radii, max_masses);
title('Scaled ADM Mass (M) vs Radial Coordinate (r)');
xlabel('r (Radius)');
ylabel('M (Scaled Mass)');
grid on;
% Plotting m vs phi
figure;
plot(phi_values, mass_values, 'color','red');
title('Scaled ADM Mass (M) vs \phi_c');
xlabel('\phi_c');
ylabel('M (Scaled Mass)');
grid on;
% Function to define the differential equations
function dydr = bsode(r, y, p)
G = 1; % Gravitational constant
omega = p(1);
m=1;
a = y(1);
alpha = y(2);
psi0 = y(3);
phi = y(4);
dydr = zeros(4, 1);
dydr(1) = (0.5) * ((a/r) * ((1 - a^2) + 4 * pi * G * r * a * ...
(psi0^2 * a^2 * (m^2 + omega^2 / alpha^2) + phi^2)));
dydr(2) = (alpha/2) * (((a^2 - 1)/r) + 4 * pi * r * ...
(psi0^2 * a^2 * (omega^2 / alpha^2 - m^2) + phi^2));
dydr(3) = phi;
dydr(4) = -(1 + a^2 - 4 * pi * r^2 * psi0^2 * a^2 * m^2) * (phi/r) - ...
(omega^2 / alpha^2 - m^2) * psi0 * a^2;
end
% Function to define boundary conditions
function res = bsbc(ya, yb, p, phi_c)
res = [ya(1) - 1; % a(0) = 1
ya(2) - 1; % alpha(0) = 1
ya(3) - phi_c; % psi0(0) = phi_c
yb(3); % psi0(infinity) = 0
ya(4)]; % phi(0) = 0
end
% Function for calculating ADM mass
function M = ADM_mass(r, a)
M = (1 - a.^-2) .* (r / 2);
end
10 Comments
T
on 12 Nov 2024
is there a different solver that i can use in this case? i want phi_c from 0.01 to 0.8. i want this plot
T
on 13 Nov 2024
Torsten
on 13 Nov 2024
Yes I am referring to a paper where this plot is given.
Could you make the plot with the code from above up to the values of phi_c where the code works ? Do the graphs agree up to this point ?
Can bvp5c be better than bvp4c?
Replace the name "bvp4c" by "bvp5c" in your code and you will see. But I doubt that the solver is the reason for your problems.
T
on 13 Nov 2024
Torsten
on 13 Nov 2024
Yes it is working till 0.15 value of phi_c and the graph is correct till here.
Can you give the command how to generate the graph up to phi_c = 0.15 ? It's different from the graphs you implemented up to now.
here!
% Define global variable for initial scalar field values
phi_c_values = logspace(log10(0.01), log10(0.15), 200); % Extended range for phi_c
% Physical constants
M_pl_squared = 1.0;
m = 1.0; % Scalar mass
% Main loop over phi_c values
omega = 0.864;
infinity = 15;
x_init = logspace(log10(1e-5), log10(infinity), 2000); % Higher resolution for x_init
% Lists to store maximum masses and radii
max_masses = [];
max_radii = [];
phi_values = [];
mass_values = [];
options = bvpset('RelTol', 1e-10, 'AbsTol', 1e-10); % Higher precision for solver
% Initial solution guess
solinit = bvpinit(x_init, [1, 1, phi_c_values(1), 0], omega);
for phi_c = phi_c_values
try
% Solve the BVP
sol = bvp4c(@bsode, @(x, y, p) bsbc(x, y, p, phi_c), solinit, options);
f = sol.y;
% Find where the solution sufficiently approaches zero for psi0
i = find(f(3, :) < 1e-5, 1);
if isempty(i)
i = length(f(3, :)); % Use last index if no value falls below threshold
end
r = sol.x(1:i);
a = f(1, 1:i);
% Scale alpha and omega
alpha_last = 1 / f(2, i);
scaled_alpha = alpha_last * f(2, 1:i);
scaled_omega = alpha_last * sol.parameters(1);
% Calculate ADM mass
M = ADM_mass(r, a);
M_max = 0.633 * M_pl_squared / m;
M_scaled = M / M_max;
% Find maximum mass and corresponding radius
max_mass = max(M_scaled);
target_mass = 0.99 * max_mass;
index = find(abs(M_scaled - target_mass) == min(abs(M_scaled - target_mass)), 1);
max_radius = r(index);
max_M = M_scaled(index);
% Store results
max_masses(end + 1) = max_M;
max_radii(end + 1) = max_radius;
phi_values(end + 1) = phi_c;
mass_values(end + 1) = max_mass;
fprintf('phi_c = %.5f, max_mass = %.5f, max_radius = %.5f\n', phi_c, max_M, max_radius);
% Update solinit for the next iteration
solinit = bvpinit(sol.x, @(x) interp1(sol.x, sol.y.', x, 'linear', 'extrap'), sol.parameters(1));
catch ME
fprintf('Solver failed for phi_c = %.5f: %s\n', phi_c, ME.message);
break; % Exit loop if solver fails at higher phi_c
end
end
% Plotting mass-radius curve
figure;
plot(max_radii, max_masses);
title('Scaled ADM Mass (M) vs Radial Coordinate (r)');
xlabel('r (Radius)');
ylabel('M (Scaled Mass)');
grid on;
% Plotting m vs phi
figure;
plot(phi_values, mass_values, 'color', 'red');
title('Scaled ADM Mass (M) vs \phi_c');
xlabel('\phi_c');
ylabel('M (Scaled Mass)');
grid on;
% Function to define the differential equations
function dydr = bsode(r, y, p)
G = 1; % Gravitational constant
omega = p(1);
m = 1;
a = y(1);
alpha = y(2);
psi0 = y(3);
phi = y(4);
dydr = zeros(4, 1);
dydr(1) = (0.5) * ((a/r) * ((1 - a^2) + 4 * pi * G * r * a * ...
(psi0^2 * a^2 * (m^2 + omega^2 / alpha^2) + phi^2)));
dydr(2) = (alpha/2) * (((a^2 - 1)/r) + 4 * pi * r * ...
(psi0^2 * a^2 * (omega^2 / alpha^2 - m^2) + phi^2));
dydr(3) = phi;
dydr(4) = -(1 + a^2 - 4 * pi * r^2 * psi0^2 * a^2 * m^2) * (phi/r) - ...
(omega^2 / alpha^2 - m^2) * psi0 * a^2;
end
% Function to define boundary conditions
function res = bsbc(ya, yb, p, phi_c)
res = [ya(1) - 1; % a(0) = 1
ya(2) - 1; % alpha(0) = 1
ya(3) - phi_c; % psi0(0) = phi_c
yb(3); % psi0(infinity) = 0
ya(4)]; % phi(0) = 0
end
% Function for calculating ADM mass
function M = ADM_mass(r, a)
M = (1 - a.^-2) .* (r / 2);
end
Torsten
on 14 Nov 2024
The graphs are different from the one you posted.
T
on 14 Nov 2024
Torsten
on 14 Nov 2024
If you get a different solution, you use different equations and/or parameters in my opinion.
Categories
Find more on Boundary Value Problems in Help Center and File Exchange
on 14 Nov 2024
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
