getting singular jacobian error while using bvp4c solver

Question

T 2024-11-12，13:32

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2165604-getting-singular-jacobian-error-while-using-bvp4c-solver

评论： Torsten 2024-11-14，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;

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

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

Answer 1

Torsten 2024-11-12，14:23

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2165604-getting-singular-jacobian-error-while-using-bvp4c-solver#answer_1544054

在 MATLAB Online 中打开

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 个评论
显示 8更早的评论隐藏 8更早的评论

T 2024-11-14，5:08

在 MATLAB Online 中打开

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

Warning: Unable to meet the tolerance without using more than 2500 mesh points.
The last mesh of 2000 points and the solution are available in the output argument.
The maximum residual is 4.21097e-07, while requested accuracy is 1e-10.

phi_c = 0.01000, max_mass = 0.05267, max_radius = 12.82686 phi_c = 0.01014, max_mass = 0.05399, max_radius = 12.82686 phi_c = 0.01028, max_mass = 0.05537, max_radius = 12.82686 phi_c = 0.01042, max_mass = 0.05679, max_radius = 12.82686 phi_c = 0.01056, max_mass = 0.05824, max_radius = 12.82686 phi_c = 0.01070, max_mass = 0.05972, max_radius = 12.82686 phi_c = 0.01085, max_mass = 0.06124, max_radius = 12.82686 phi_c = 0.01100, max_mass = 0.06280, max_radius = 12.82686 phi_c = 0.01115, max_mass = 0.06439, max_radius = 12.82686 phi_c = 0.01130, max_mass = 0.06601, max_radius = 12.82686 phi_c = 0.01146, max_mass = 0.06768, max_radius = 12.82686 phi_c = 0.01161, max_mass = 0.06938, max_radius = 12.82686 phi_c = 0.01177, max_mass = 0.07112, max_radius = 12.82686 phi_c = 0.01194, max_mass = 0.07290, max_radius = 12.82686 phi_c = 0.01210, max_mass = 0.07471, max_radius = 12.82686 phi_c = 0.01226, max_mass = 0.07657, max_radius = 12.82686 phi_c = 0.01243, max_mass = 0.07847, max_radius = 12.82686 phi_c = 0.01260, max_mass = 0.08041, max_radius = 12.82686 phi_c = 0.01278, max_mass = 0.08240, max_radius = 12.82686 phi_c = 0.01295, max_mass = 0.08442, max_radius = 12.82686 phi_c = 0.01313, max_mass = 0.08650, max_radius = 12.82686 phi_c = 0.01331, max_mass = 0.08861, max_radius = 12.82686 phi_c = 0.01349, max_mass = 0.09077, max_radius = 12.82686 phi_c = 0.01368, max_mass = 0.09297, max_radius = 12.82686 phi_c = 0.01386, max_mass = 0.09522, max_radius = 12.82686 phi_c = 0.01405, max_mass = 0.09752, max_radius = 12.82686 phi_c = 0.01424, max_mass = 0.09976, max_radius = 12.73593 phi_c = 0.01444, max_mass = 0.10215, max_radius = 12.73593 phi_c = 0.01464, max_mass = 0.10459, max_radius = 12.73593 phi_c = 0.01484, max_mass = 0.10708, max_radius = 12.73593 phi_c = 0.01504, max_mass = 0.10961, max_radius = 12.73593 phi_c = 0.01525, max_mass = 0.11220, max_radius = 12.73593 phi_c = 0.01546, max_mass = 0.11484, max_radius = 12.73593 phi_c = 0.01567, max_mass = 0.11752, max_radius = 12.73593 phi_c = 0.01588, max_mass = 0.12026, max_radius = 12.73593 phi_c = 0.01610, max_mass = 0.12305, max_radius = 12.73593 phi_c = 0.01632, max_mass = 0.12589, max_radius = 12.73593 phi_c = 0.01655, max_mass = 0.12879, max_radius = 12.73593 phi_c = 0.01677, max_mass = 0.13174, max_radius = 12.73593 phi_c = 0.01700, max_mass = 0.13474, max_radius = 12.73593 phi_c = 0.01723, max_mass = 0.13779, max_radius = 12.73593 phi_c = 0.01747, max_mass = 0.14090, max_radius = 12.73593 phi_c = 0.01771, max_mass = 0.14391, max_radius = 12.64565 phi_c = 0.01795, max_mass = 0.14712, max_radius = 12.64565 phi_c = 0.01820, max_mass = 0.15039, max_radius = 12.64565 phi_c = 0.01845, max_mass = 0.15371, max_radius = 12.64565 phi_c = 0.01870, max_mass = 0.15708, max_radius = 12.64565 phi_c = 0.01896, max_mass = 0.16051, max_radius = 12.64565 phi_c = 0.01922, max_mass = 0.16400, max_radius = 12.64565 phi_c = 0.01948, max_mass = 0.16753, max_radius = 12.64565 phi_c = 0.01975, max_mass = 0.17112, max_radius = 12.64565 phi_c = 0.02002, max_mass = 0.17477, max_radius = 12.64565 phi_c = 0.02029, max_mass = 0.17847, max_radius = 12.64565 phi_c = 0.02057, max_mass = 0.18204, max_radius = 12.55601 phi_c = 0.02085, max_mass = 0.18584, max_radius = 12.55601 phi_c = 0.02114, max_mass = 0.18970, max_radius = 12.55601 phi_c = 0.02143, max_mass = 0.19361, max_radius = 12.55601 phi_c = 0.02172, max_mass = 0.19757, max_radius = 12.55601 phi_c = 0.02202, max_mass = 0.20158, max_radius = 12.55601 phi_c = 0.02232, max_mass = 0.20564, max_radius = 12.55601 phi_c = 0.02263, max_mass = 0.20975, max_radius = 12.55601 phi_c = 0.02294, max_mass = 0.21370, max_radius = 12.46700 phi_c = 0.02325, max_mass = 0.21791, max_radius = 12.46700 phi_c = 0.02357, max_mass = 0.22217, max_radius = 12.46700 phi_c = 0.02389, max_mass = 0.22647, max_radius = 12.46700 phi_c = 0.02422, max_mass = 0.23082, max_radius = 12.46700 phi_c = 0.02455, max_mass = 0.23521, max_radius = 12.46700 phi_c = 0.02489, max_mass = 0.23965, max_radius = 12.46700 phi_c = 0.02523, max_mass = 0.24390, max_radius = 12.37862 phi_c = 0.02557, max_mass = 0.24842, max_radius = 12.37862 phi_c = 0.02592, max_mass = 0.25299, max_radius = 12.37862 phi_c = 0.02628, max_mass = 0.25759, max_radius = 12.37862 phi_c = 0.02664, max_mass = 0.26222, max_radius = 12.37862 phi_c = 0.02700, max_mass = 0.26665, max_radius = 12.29087 phi_c = 0.02737, max_mass = 0.27136, max_radius = 12.29087 phi_c = 0.02775, max_mass = 0.27611, max_radius = 12.29087 phi_c = 0.02813, max_mass = 0.28088, max_radius = 12.29087 phi_c = 0.02851, max_mass = 0.28568, max_radius = 12.29087 phi_c = 0.02891, max_mass = 0.29026, max_radius = 12.20375 phi_c = 0.02930, max_mass = 0.29512, max_radius = 12.20375 phi_c = 0.02970, max_mass = 0.30000, max_radius = 12.20375 phi_c = 0.03011, max_mass = 0.30491, max_radius = 12.20375 phi_c = 0.03052, max_mass = 0.30984, max_radius = 12.20375 phi_c = 0.03094, max_mass = 0.31452, max_radius = 12.11724 phi_c = 0.03136, max_mass = 0.31949, max_radius = 12.11724 phi_c = 0.03179, max_mass = 0.32448, max_radius = 12.11724 phi_c = 0.03223, max_mass = 0.32919, max_radius = 12.03134 phi_c = 0.03267, max_mass = 0.33421, max_radius = 12.03134 phi_c = 0.03312, max_mass = 0.33923, max_radius = 12.03134 phi_c = 0.03357, max_mass = 0.34426, max_radius = 12.03134 phi_c = 0.03403, max_mass = 0.34901, max_radius = 11.94605 phi_c = 0.03450, max_mass = 0.35405, max_radius = 11.94605 phi_c = 0.03497, max_mass = 0.35909, max_radius = 11.94605 phi_c = 0.03545, max_mass = 0.36384, max_radius = 11.86137 phi_c = 0.03594, max_mass = 0.36888, max_radius = 11.86137 phi_c = 0.03643, max_mass = 0.37392, max_radius = 11.86137 phi_c = 0.03693, max_mass = 0.37866, max_radius = 11.77728 phi_c = 0.03743, max_mass = 0.38369, max_radius = 11.77728 phi_c = 0.03795, max_mass = 0.38870, max_radius = 11.77728 phi_c = 0.03847, max_mass = 0.39341, max_radius = 11.69380 phi_c = 0.03899, max_mass = 0.39840, max_radius = 11.69380 phi_c = 0.03953, max_mass = 0.40308, max_radius = 11.61090 phi_c = 0.04007, max_mass = 0.40805, max_radius = 11.61090 phi_c = 0.04062, max_mass = 0.41268, max_radius = 11.52860 phi_c = 0.04118, max_mass = 0.41761, max_radius = 11.52860 phi_c = 0.04174, max_mass = 0.42221, max_radius = 11.44687 phi_c = 0.04231, max_mass = 0.42709, max_radius = 11.44687 phi_c = 0.04289, max_mass = 0.43195, max_radius = 11.44687 phi_c = 0.04348, max_mass = 0.43647, max_radius = 11.36573 phi_c = 0.04407, max_mass = 0.44097, max_radius = 11.28516 phi_c = 0.04468, max_mass = 0.44574, max_radius = 11.28516 phi_c = 0.04529, max_mass = 0.45018, max_radius = 11.20516 phi_c = 0.04591, max_mass = 0.45489, max_radius = 11.20516 phi_c = 0.04654, max_mass = 0.45927, max_radius = 11.12573 phi_c = 0.04718, max_mass = 0.46392, max_radius = 11.12573 phi_c = 0.04782, max_mass = 0.46822, max_radius = 11.04686 phi_c = 0.04848, max_mass = 0.47249, max_radius = 10.96855 phi_c = 0.04914, max_mass = 0.47704, max_radius = 10.96855 phi_c = 0.04982, max_mass = 0.48123, max_radius = 10.89080 phi_c = 0.05050, max_mass = 0.48540, max_radius = 10.81359 phi_c = 0.05119, max_mass = 0.48952, max_radius = 10.73694 phi_c = 0.05189, max_mass = 0.49390, max_radius = 10.73694 phi_c = 0.05260, max_mass = 0.49795, max_radius = 10.66083 phi_c = 0.05332, max_mass = 0.50195, max_radius = 10.58526 phi_c = 0.05406, max_mass = 0.50592, max_radius = 10.51022 phi_c = 0.05480, max_mass = 0.51014, max_radius = 10.51022 phi_c = 0.05555, max_mass = 0.51402, max_radius = 10.43571 phi_c = 0.05631, max_mass = 0.51786, max_radius = 10.36174 phi_c = 0.05708, max_mass = 0.52165, max_radius = 10.28828 phi_c = 0.05786, max_mass = 0.52541, max_radius = 10.21535 phi_c = 0.05865, max_mass = 0.52912, max_radius = 10.14294 phi_c = 0.05946, max_mass = 0.53279, max_radius = 10.07104 phi_c = 0.06027, max_mass = 0.53641, max_radius = 9.99965 phi_c = 0.06110, max_mass = 0.53999, max_radius = 9.92876 phi_c = 0.06194, max_mass = 0.54352, max_radius = 9.85838 phi_c = 0.06278, max_mass = 0.54701, max_radius = 9.78849 phi_c = 0.06364, max_mass = 0.55046, max_radius = 9.71910 phi_c = 0.06452, max_mass = 0.55357, max_radius = 9.58180 phi_c = 0.06540, max_mass = 0.55693, max_radius = 9.51388 phi_c = 0.06630, max_mass = 0.56024, max_radius = 9.44643 phi_c = 0.06720, max_mass = 0.56350, max_radius = 9.37947 phi_c = 0.06813, max_mass = 0.56671, max_radius = 9.31298 phi_c = 0.06906, max_mass = 0.56959, max_radius = 9.18141 phi_c = 0.07001, max_mass = 0.57271, max_radius = 9.11633 phi_c = 0.07096, max_mass = 0.57578, max_radius = 9.05170 phi_c = 0.07194, max_mass = 0.57851, max_radius = 8.92383 phi_c = 0.07292, max_mass = 0.58148, max_radius = 8.86057 phi_c = 0.07392, max_mass = 0.58411, max_radius = 8.73539 phi_c = 0.07493, max_mass = 0.58697, max_radius = 8.67347 phi_c = 0.07596, max_mass = 0.58957, max_radius = 8.56631 phi_c = 0.07700, max_mass = 0.59218, max_radius = 8.47538 phi_c = 0.07806, max_mass = 0.59474, max_radius = 8.38542 phi_c = 0.07913, max_mass = 0.59723, max_radius = 8.29641 phi_c = 0.08021, max_mass = 0.59967, max_radius = 8.20835 phi_c = 0.08131, max_mass = 0.60204, max_radius = 8.12123 phi_c = 0.08242, max_mass = 0.60434, max_radius = 8.03503 phi_c = 0.08355, max_mass = 0.60658, max_radius = 7.94974 phi_c = 0.08470, max_mass = 0.60875, max_radius = 7.86536 phi_c = 0.08586, max_mass = 0.61065, max_radius = 7.74033 phi_c = 0.08703, max_mass = 0.61267, max_radius = 7.65822 phi_c = 0.08823, max_mass = 0.61461, max_radius = 7.57688 phi_c = 0.08944, max_mass = 0.61628, max_radius = 7.45670 phi_c = 0.09066, max_mass = 0.61806, max_radius = 7.37755 phi_c = 0.09190, max_mass = 0.61975, max_radius = 7.29924 phi_c = 0.09316, max_mass = 0.62115, max_radius = 7.18326 phi_c = 0.09444, max_mass = 0.62265, max_radius = 7.10697 phi_c = 0.09573, max_mass = 0.62385, max_radius = 6.99424 phi_c = 0.09704, max_mass = 0.62515, max_radius = 6.92000 phi_c = 0.09837, max_mass = 0.62613, max_radius = 6.80999 phi_c = 0.09972, max_mass = 0.62720, max_radius = 6.73775 phi_c = 0.10109, max_mass = 0.62795, max_radius = 6.63084 phi_c = 0.10247, max_mass = 0.62876, max_radius = 6.56045 phi_c = 0.10388, max_mass = 0.62924, max_radius = 6.45621 phi_c = 0.10530, max_mass = 0.62958, max_radius = 6.35376 phi_c = 0.10674, max_mass = 0.62996, max_radius = 6.28632 phi_c = 0.10821, max_mass = 0.62999, max_radius = 6.18643 phi_c = 0.10969, max_mass = 0.62984, max_radius = 6.08826 phi_c = 0.11119, max_mass = 0.62952, max_radius = 5.99148 phi_c = 0.11271, max_mass = 0.62900, max_radius = 5.89644 phi_c = 0.11426, max_mass = 0.62839, max_radius = 5.82339 phi_c = 0.11582, max_mass = 0.62749, max_radius = 5.73867 phi_c = 0.11741, max_mass = 0.62631, max_radius = 5.64758 phi_c = 0.11902, max_mass = 0.62488, max_radius = 5.56027 phi_c = 0.12065, max_mass = 0.62316, max_radius = 5.47204

Warning: Unable to meet the tolerance without using more than 2500 mesh points.
The last mesh of 1908 points and the solution are available in the output argument.
The maximum residual is 1.70477e-10, while requested accuracy is 1e-10.

phi_c = 0.12230, max_mass = 0.62114, max_radius = 5.38506

Warning: Unable to meet the tolerance without using more than 2500 mesh points.
The last mesh of 1908 points and the solution are available in the output argument.
The maximum residual is 2.70565e-10, while requested accuracy is 1e-10.

phi_c = 0.12398, max_mass = 0.61881, max_radius = 5.30202

Warning: Unable to meet the tolerance without using more than 2500 mesh points.
The last mesh of 1908 points and the solution are available in the output argument.
The maximum residual is 4.32897e-10, while requested accuracy is 1e-10.

phi_c = 0.12568, max_mass = 0.61597, max_radius = 5.19688

Warning: Unable to meet the tolerance without using more than 2500 mesh points.
The last mesh of 1908 points and the solution are available in the output argument.
The maximum residual is 7.02202e-10, while requested accuracy is 1e-10.

phi_c = 0.12740, max_mass = 0.61289, max_radius = 5.11674

Warning: Unable to meet the tolerance without using more than 2500 mesh points.
The last mesh of 1908 points and the solution are available in the output argument.
The maximum residual is 1.15997e-09, while requested accuracy is 1e-10.

phi_c = 0.12915, max_mass = 0.60935, max_radius = 5.03546

Warning: Unable to meet the tolerance without using more than 2500 mesh points.
The last mesh of 1908 points and the solution are available in the output argument.
The maximum residual is 1.96208e-09, while requested accuracy is 1e-10.

phi_c = 0.13092, max_mass = 0.60534, max_radius = 4.95769

Warning: Unable to meet the tolerance without using more than 2500 mesh points.
The last mesh of 1908 points and the solution are available in the output argument.
The maximum residual is 3.41978e-09, while requested accuracy is 1e-10.

phi_c = 0.13271, max_mass = 0.60064, max_radius = 4.85949

Warning: Unable to meet the tolerance without using more than 2500 mesh points.
The last mesh of 1908 points and the solution are available in the output argument.
The maximum residual is 6.19073e-09, while requested accuracy is 1e-10.

phi_c = 0.13453, max_mass = 0.59541, max_radius = 4.78021

Warning: Unable to meet the tolerance without using more than 2500 mesh points.
The last mesh of 1908 points and the solution are available in the output argument.
The maximum residual is 1.17598e-08, while requested accuracy is 1e-10.

phi_c = 0.13637, max_mass = 0.58947, max_radius = 4.70422

Warning: Unable to meet the tolerance without using more than 2500 mesh points.
The last mesh of 1908 points and the solution are available in the output argument.
The maximum residual is 2.37693e-08, while requested accuracy is 1e-10.

phi_c = 0.13824, max_mass = 0.58252, max_radius = 4.60906

Warning: Unable to meet the tolerance without using more than 2500 mesh points.
The last mesh of 1908 points and the solution are available in the output argument.
The maximum residual is 5.21514e-08, while requested accuracy is 1e-10.

phi_c = 0.14013, max_mass = 0.57468, max_radius = 4.53983

Warning: Unable to meet the tolerance without using more than 2500 mesh points.
The last mesh of 1908 points and the solution are available in the output argument.
The maximum residual is 1.28017e-07, while requested accuracy is 1e-10.

phi_c = 0.14205, max_mass = 0.56545, max_radius = 4.46186 phi_c = 0.14400, max_mass = 0.55450, max_radius = 4.38907

Warning: Unable to meet the tolerance without using more than 2500 mesh points.
The last mesh of 1810 points and the solution are available in the output argument.
The maximum residual is 5.76537e-10, while requested accuracy is 1e-10.

phi_c = 0.14597, max_mass = 0.54116, max_radius = 4.33668

Warning: Unable to meet the tolerance without using more than 2500 mesh points.
The last mesh of 1810 points and the solution are available in the output argument.
The maximum residual is 1.42595e-08, while requested accuracy is 1e-10.

phi_c = 0.14797, max_mass = 0.52372, max_radius = 4.28868

Warning: Unable to meet the tolerance without using more than 2500 mesh points.
The last mesh of 1810 points and the solution are available in the output argument.
The maximum residual is 1.03788e-07, while requested accuracy is 1e-10.

phi_c = 0.15000, max_mass = 0.49679, max_radius = 4.33668

% 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

T 2024-11-14，12:20

yes a bit different. but this is the best i can get from this solver

Torsten 2024-11-14，12:30

If you get a different solution, you use different equations and/or parameters in my opinion.

请先登录，再进行评论。

getting singular jacobian error while using bvp4c solver

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

回答（1 个）

10 个评论
显示 8更早的评论隐藏 8更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

getting singular jacobian error while using bvp4c solver

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

回答（1 个）

10 个评论 显示 8更早的评论隐藏 8更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

10 个评论
显示 8更早的评论隐藏 8更早的评论