getting singular jacobian error while using bvp4c solver

40 views (last 30 days)
T
T on 12 Nov 2024 at 13:32
Commented: Torsten on 14 Nov 2024 at 12:30
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)

Torsten
Torsten on 12 Nov 2024 at 14:23
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
phi_c = 0.01000, max_mass = 0.05268, max_radius = 12.89409 phi_c = 0.01070, max_mass = 0.05977, max_radius = 12.89409 phi_c = 0.01141, max_mass = 0.06718, max_radius = 12.89409 phi_c = 0.01211, max_mass = 0.07490, max_radius = 12.89409 phi_c = 0.01281, max_mass = 0.08291, max_radius = 12.89409 phi_c = 0.01352, max_mass = 0.09116, max_radius = 12.89409 phi_c = 0.01422, max_mass = 0.09965, max_radius = 12.89409 phi_c = 0.01492, max_mass = 0.10817, max_radius = 12.74978 phi_c = 0.01563, max_mass = 0.11703, max_radius = 12.74978 phi_c = 0.01633, max_mass = 0.12605, max_radius = 12.74978 phi_c = 0.01704, max_mass = 0.13520, max_radius = 12.74978 phi_c = 0.01774, max_mass = 0.14446, max_radius = 12.74978 phi_c = 0.01844, max_mass = 0.15381, max_radius = 12.74978 phi_c = 0.01915, max_mass = 0.16323, max_radius = 12.74978 phi_c = 0.01985, max_mass = 0.17269, max_radius = 12.74978 phi_c = 0.02055, max_mass = 0.18218, max_radius = 12.74978 phi_c = 0.02126, max_mass = 0.19169, max_radius = 12.74978 phi_c = 0.02196, max_mass = 0.20118, max_radius = 12.74978 phi_c = 0.02266, max_mass = 0.20967, max_radius = 12.31686 phi_c = 0.02337, max_mass = 0.21909, max_radius = 12.31686 phi_c = 0.02407, max_mass = 0.22846, max_radius = 12.31686 phi_c = 0.02477, max_mass = 0.23777, max_radius = 12.31686 phi_c = 0.02548, max_mass = 0.24700, max_radius = 12.31686 phi_c = 0.02618, max_mass = 0.25615, max_radius = 12.31686 phi_c = 0.02688, max_mass = 0.26520, max_radius = 12.31686 phi_c = 0.02759, max_mass = 0.27415, max_radius = 12.31686 phi_c = 0.02829, max_mass = 0.28298, max_radius = 12.31686 phi_c = 0.02899, max_mass = 0.29169, max_radius = 12.31686 phi_c = 0.02970, max_mass = 0.30028, max_radius = 12.31686 phi_c = 0.03040, max_mass = 0.30874, max_radius = 12.31686 phi_c = 0.03111, max_mass = 0.31706, max_radius = 12.31686 phi_c = 0.03181, max_mass = 0.32523, max_radius = 12.31686 phi_c = 0.03251, max_mass = 0.33189, max_radius = 11.88394 phi_c = 0.03322, max_mass = 0.33979, max_radius = 11.88394 phi_c = 0.03392, max_mass = 0.34755, max_radius = 11.88394 phi_c = 0.03462, max_mass = 0.35516, max_radius = 11.88394 phi_c = 0.03533, max_mass = 0.36262, max_radius = 11.88394 phi_c = 0.03603, max_mass = 0.36992, max_radius = 11.88394 phi_c = 0.03673, max_mass = 0.37708, max_radius = 11.88394 phi_c = 0.03744, max_mass = 0.38331, max_radius = 11.66748 phi_c = 0.03814, max_mass = 0.39018, max_radius = 11.66748 phi_c = 0.03884, max_mass = 0.39690, max_radius = 11.66748 phi_c = 0.03955, max_mass = 0.40347, max_radius = 11.66748 phi_c = 0.04025, max_mass = 0.40989, max_radius = 11.66748 phi_c = 0.04095, max_mass = 0.41537, max_radius = 11.45102 phi_c = 0.04166, max_mass = 0.42152, max_radius = 11.45102 phi_c = 0.04236, max_mass = 0.42753, max_radius = 11.45102 phi_c = 0.04307, max_mass = 0.43340, max_radius = 11.45102 phi_c = 0.04377, max_mass = 0.43913, max_radius = 11.45102 phi_c = 0.04447, max_mass = 0.44473, max_radius = 11.45102 phi_c = 0.04518, max_mass = 0.44898, max_radius = 11.12633 phi_c = 0.04588, max_mass = 0.45434, max_radius = 11.12633 phi_c = 0.04658, max_mass = 0.45958, max_radius = 11.12633 phi_c = 0.04729, max_mass = 0.46470, max_radius = 11.12633 phi_c = 0.04799, max_mass = 0.46969, max_radius = 11.12633 phi_c = 0.04869, max_mass = 0.47456, max_radius = 11.12633 phi_c = 0.04940, max_mass = 0.47807, max_radius = 10.80164 phi_c = 0.05010, max_mass = 0.48274, max_radius = 10.80164 phi_c = 0.05080, max_mass = 0.48731, max_radius = 10.80164 phi_c = 0.05151, max_mass = 0.49176, max_radius = 10.80164 phi_c = 0.05221, max_mass = 0.49610, max_radius = 10.80164 phi_c = 0.05291, max_mass = 0.50034, max_radius = 10.80164 phi_c = 0.05362, max_mass = 0.50323, max_radius = 10.47695 phi_c = 0.05432, max_mass = 0.50731, max_radius = 10.47695 phi_c = 0.05503, max_mass = 0.51129, max_radius = 10.47695 phi_c = 0.05573, max_mass = 0.51452, max_radius = 10.31460 phi_c = 0.05643, max_mass = 0.51833, max_radius = 10.31460 phi_c = 0.05714, max_mass = 0.52205, max_radius = 10.31460 phi_c = 0.05784, max_mass = 0.52503, max_radius = 10.15226 phi_c = 0.05854, max_mass = 0.52860, max_radius = 10.15226 phi_c = 0.05925, max_mass = 0.53208, max_radius = 10.15226 phi_c = 0.05995, max_mass = 0.53448, max_radius = 9.90874 phi_c = 0.06065, max_mass = 0.53783, max_radius = 9.90874 phi_c = 0.06136, max_mass = 0.54110, max_radius = 9.90874 phi_c = 0.06206, max_mass = 0.54429, max_radius = 9.90874 phi_c = 0.06276, max_mass = 0.54640, max_radius = 9.66522 phi_c = 0.06347, max_mass = 0.54947, max_radius = 9.66522 phi_c = 0.06417, max_mass = 0.55247, max_radius = 9.66522 phi_c = 0.06487, max_mass = 0.55540, max_radius = 9.66522 phi_c = 0.06558, max_mass = 0.55725, max_radius = 9.42170 phi_c = 0.06628, max_mass = 0.56007, max_radius = 9.42170 phi_c = 0.06698, max_mass = 0.56282, max_radius = 9.42170 phi_c = 0.06769, max_mass = 0.56551, max_radius = 9.42170 phi_c = 0.06839, max_mass = 0.56712, max_radius = 9.17819 phi_c = 0.06910, max_mass = 0.56971, max_radius = 9.17819 phi_c = 0.06980, max_mass = 0.57224, max_radius = 9.17819 phi_c = 0.07050, max_mass = 0.57471, max_radius = 9.17819 phi_c = 0.07121, max_mass = 0.57609, max_radius = 8.93467 phi_c = 0.07191, max_mass = 0.57847, max_radius = 8.93467 phi_c = 0.07261, max_mass = 0.58079, max_radius = 8.93467 phi_c = 0.07332, max_mass = 0.58254, max_radius = 8.81291 phi_c = 0.07402, max_mass = 0.58476, max_radius = 8.81291 phi_c = 0.07472, max_mass = 0.58613, max_radius = 8.63027 phi_c = 0.07543, max_mass = 0.58827, max_radius = 8.63027 phi_c = 0.07613, max_mass = 0.59036, max_radius = 8.63027 phi_c = 0.07683, max_mass = 0.59157, max_radius = 8.44763 phi_c = 0.07754, max_mass = 0.59358, max_radius = 8.44763 phi_c = 0.07824, max_mass = 0.59552, max_radius = 8.44763 phi_c = 0.07894, max_mass = 0.59660, max_radius = 8.26500 phi_c = 0.07965, max_mass = 0.59847, max_radius = 8.26500 phi_c = 0.08035, max_mass = 0.60029, max_radius = 8.26500 phi_c = 0.08106, max_mass = 0.60121, max_radius = 8.08236 phi_c = 0.08176, max_mass = 0.60296, max_radius = 8.08236 phi_c = 0.08246, max_mass = 0.60465, max_radius = 8.08236 phi_c = 0.08317, max_mass = 0.60544, max_radius = 7.89972 phi_c = 0.08387, max_mass = 0.60706, max_radius = 7.89972 phi_c = 0.08457, max_mass = 0.60864, max_radius = 7.89972 phi_c = 0.08528, max_mass = 0.60928, max_radius = 7.71708 phi_c = 0.08598, max_mass = 0.61079, max_radius = 7.71708 phi_c = 0.08668, max_mass = 0.61224, max_radius = 7.71708 phi_c = 0.08739, max_mass = 0.61365, max_radius = 7.71708 phi_c = 0.08809, max_mass = 0.61414, max_radius = 7.53444 phi_c = 0.08879, max_mass = 0.61548, max_radius = 7.53444 phi_c = 0.08950, max_mass = 0.61655, max_radius = 7.48878 phi_c = 0.09020, max_mass = 0.61711, max_radius = 7.35180 phi_c = 0.09090, max_mass = 0.61834, max_radius = 7.35180 phi_c = 0.09161, max_mass = 0.61953, max_radius = 7.35180 phi_c = 0.09231, max_mass = 0.61997, max_radius = 7.21482 phi_c = 0.09302, max_mass = 0.62108, max_radius = 7.21482 phi_c = 0.09372, max_mass = 0.62214, max_radius = 7.21482 phi_c = 0.09442, max_mass = 0.62247, max_radius = 7.07785 phi_c = 0.09513, max_mass = 0.62346, max_radius = 7.07785 phi_c = 0.09583, max_mass = 0.62369, max_radius = 6.94087 phi_c = 0.09653, max_mass = 0.62461, max_radius = 6.94087 phi_c = 0.09724, max_mass = 0.62549, max_radius = 6.94087 phi_c = 0.09794, max_mass = 0.62559, max_radius = 6.80389 phi_c = 0.09864, max_mass = 0.62640, max_radius = 6.80389 phi_c = 0.09935, max_mass = 0.62716, max_radius = 6.80389 phi_c = 0.10005, max_mass = 0.62714, max_radius = 6.66691 phi_c = 0.10075, max_mass = 0.62783, max_radius = 6.66691 phi_c = 0.10146, max_mass = 0.62848, max_radius = 6.66691 phi_c = 0.10216, max_mass = 0.62833, max_radius = 6.52993 phi_c = 0.10286, max_mass = 0.62890, max_radius = 6.52993 phi_c = 0.10357, max_mass = 0.62943, max_radius = 6.52993 phi_c = 0.10427, max_mass = 0.62916, max_radius = 6.39295 phi_c = 0.10497, max_mass = 0.62961, max_radius = 6.39295 phi_c = 0.10568, max_mass = 0.62963, max_radius = 6.32446 phi_c = 0.10638, max_mass = 0.63000, max_radius = 6.32446 phi_c = 0.10709, max_mass = 0.62973, max_radius = 6.21806 phi_c = 0.10779, max_mass = 0.63002, max_radius = 6.21806 phi_c = 0.10849, max_mass = 0.62964, max_radius = 6.11167 phi_c = 0.10920, max_mass = 0.62986, max_radius = 6.11167 phi_c = 0.10990, max_mass = 0.63002, max_radius = 6.11167 phi_c = 0.11060, max_mass = 0.62951, max_radius = 6.00344 phi_c = 0.11131, max_mass = 0.62960, max_radius = 6.00344 phi_c = 0.11201, max_mass = 0.62915, max_radius = 5.92226 phi_c = 0.11271, max_mass = 0.62915, max_radius = 5.92226 phi_c = 0.11342, max_mass = 0.62862, max_radius = 5.84109 phi_c = 0.11412, max_mass = 0.62852, max_radius = 5.84109 phi_c = 0.11482, max_mass = 0.62783, max_radius = 5.74977 phi_c = 0.11553, max_mass = 0.62765, max_radius = 5.74977 phi_c = 0.11623, max_mass = 0.62742, max_radius = 5.74977 phi_c = 0.11693, max_mass = 0.62659, max_radius = 5.65845 phi_c = 0.11764, max_mass = 0.62626, max_radius = 5.65845 phi_c = 0.11834, max_mass = 0.62533, max_radius = 5.56713 phi_c = 0.11905, max_mass = 0.62490, max_radius = 5.56713 phi_c = 0.11975, max_mass = 0.62386, max_radius = 5.47582 phi_c = 0.12045, max_mass = 0.62334, max_radius = 5.47582 phi_c = 0.12116, max_mass = 0.62276, max_radius = 5.47582 phi_c = 0.12186, max_mass = 0.62156, max_radius = 5.38450 phi_c = 0.12256, max_mass = 0.62087, max_radius = 5.38450 phi_c = 0.12327, max_mass = 0.61956, max_radius = 5.29318 phi_c = 0.12397, max_mass = 0.61876, max_radius = 5.29318 phi_c = 0.12467, max_mass = 0.61790, max_radius = 5.29318 phi_c = 0.12538, max_mass = 0.61641, max_radius = 5.20186 phi_c = 0.12608, max_mass = 0.61543, max_radius = 5.20186 phi_c = 0.12678, max_mass = 0.61381, max_radius = 5.11054 phi_c = 0.12749, max_mass = 0.61271, max_radius = 5.11054 phi_c = 0.12819, max_mass = 0.61124, max_radius = 5.06488 phi_c = 0.12889, max_mass = 0.61000, max_radius = 5.06488 phi_c = 0.12960, max_mass = 0.60825, max_radius = 4.99639 phi_c = 0.13030, max_mass = 0.60686, max_radius = 4.99639 phi_c = 0.13101, max_mass = 0.60496, max_radius = 4.92790 phi_c = 0.13171, max_mass = 0.60341, max_radius = 4.92790 phi_c = 0.13241, max_mass = 0.60135, max_radius = 4.85941 phi_c = 0.13312, max_mass = 0.59963, max_radius = 4.85941 phi_c = 0.13382, max_mass = 0.59739, max_radius = 4.79092 phi_c = 0.13452, max_mass = 0.59549, max_radius = 4.79092 phi_c = 0.13523, max_mass = 0.59306, max_radius = 4.72243 phi_c = 0.13593, max_mass = 0.59096, max_radius = 4.72243 phi_c = 0.13663, max_mass = 0.58874, max_radius = 4.72243 phi_c = 0.13734, max_mass = 0.58598, max_radius = 4.65394 phi_c = 0.13804, max_mass = 0.58352, max_radius = 4.65394 phi_c = 0.13874, max_mass = 0.58051, max_radius = 4.58545 phi_c = 0.13945, max_mass = 0.57777, max_radius = 4.58545 phi_c = 0.14015, max_mass = 0.57446, max_radius = 4.51696 phi_c = 0.14085, max_mass = 0.57141, max_radius = 4.51696 phi_c = 0.14156, max_mass = 0.56818, max_radius = 4.51696 phi_c = 0.14226, max_mass = 0.56433, max_radius = 4.44847 phi_c = 0.14296, max_mass = 0.56068, max_radius = 4.44847 phi_c = 0.14367, max_mass = 0.55637, max_radius = 4.37999 phi_c = 0.14437, max_mass = 0.55220, max_radius = 4.37999 phi_c = 0.14508, max_mass = 0.54770, max_radius = 4.37999 phi_c = 0.14578, max_mass = 0.54250, max_radius = 4.32362 phi_c = 0.14648, max_mass = 0.53718, max_radius = 4.32362 phi_c = 0.14719, max_mass = 0.53131, max_radius = 4.32362 phi_c = 0.14789, max_mass = 0.52443, max_radius = 4.26725 phi_c = 0.14859, max_mass = 0.51691, max_radius = 4.26725 phi_c = 0.14930, max_mass = 0.50826, max_radius = 4.32362 phi_c = 0.15000, max_mass = 0.49673, max_radius = 4.32362
% 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
T on 14 Nov 2024 at 12:20
yes a bit different. but this is the best i can get from this solver
Torsten
Torsten on 14 Nov 2024 at 12:30
If you get a different solution, you use different equations and/or parameters in my opinion.

Sign in to comment.

Categories

Find more on Fluid Dynamics in Help Center and File Exchange

Products


Release

R2024a

Community Treasure Hunt

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

Start Hunting!