getting singular jacobian error while using bvp4c solver
Ältere Kommentare anzeigen
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;
Antworten (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 Kommentare
T
am 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
am 13 Nov. 2024
Torsten
am 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
am 13 Nov. 2024
Torsten
am 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
am 14 Nov. 2024
The graphs are different from the one you posted.
T
am 14 Nov. 2024
Torsten
am 14 Nov. 2024
If you get a different solution, you use different equations and/or parameters in my opinion.
Kategorien
Mehr zu Boundary Value Problems finden Sie in Hilfe-Center und File Exchange
am 14 Nov. 2024
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
