Solving bvp using shooting method to solve coupled first order odes

Question

T 2024-9-3

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2149754-solving-bvp-using-shooting-method-to-solve-coupled-first-order-odes

编辑： T 2024-9-18

I was initially trying using python but cant solve using matlab. please help me out. I have attached a sample code and the equations and boundary conditions i need to use to plot mass-radius plot for one star.

scaled equations: dadr1 = (1/2) * a * ((1 - a**2)/r1 + 4 * np.pi * G * r1 *

(psi0_1**2 * a**2 * (1 + 1 / alpha1**2)/(4*np.pi*G) + phi**2))

dalphadr1 = (alpha1/2) * (((a**2 - 1)/r1 + 4 * np.pi * G * r1 *

(psi0_1**2 * a**2 * (1 / alpha1**2 - 1)/(4*np.pi*G) + phi**2)))

dpsidr1 = phi

dphidr1 = -(1 + a**2 - r1**2 * psi0_1**2 * a**2) * (phi/r1) - \

(1 / alpha1**2 - 1)/(4*np.pi*G)**0.5 * psi0_1 * a**2

boundary conditions:ya[0] - 1, # a(0) = 1

ya[1] - 1/omega, # alpha(0) = 1 scaled

ya[2] - phi_c/(4*np.pi*G)**(1/2), # psi0(0) = phi_c scaled

yb[2] - 0, # Allow psi0(inf) to be close to zero, but not exact

yb[3] - 0

onestar2()

Warning: Unable to meet the tolerance without using more than 2500 mesh points.
The last mesh of 1500 points and the solution are available in the output argument.
The maximum error is 8.25648e-05, while requested accuracy is 1e-06.

The solution was obtained on a mesh of 1500 points. The maximum error is 8.256e-05. There were 77959 calls to the ODE function. There were 60 calls to the BC function. The parameter is 1.214 The final omega value is 0.899

Warning: Ignoring extra legend entries.

Warning: Unable to meet the tolerance without using more than 2500 mesh points.
The last mesh of 1500 points and the solution are available in the output argument.
The maximum error is 2.2136e-05, while requested accuracy is 1e-06.

The solution was obtained on a mesh of 1500 points. The maximum error is 2.214e-05. There were 78134 calls to the ODE function. There were 60 calls to the BC function. The parameter is 1.106 The final omega value is 0.948

Warning: Ignoring extra legend entries.

Warning: Unable to meet the tolerance without using more than 2500 mesh points.
The last mesh of 1500 points and the solution are available in the output argument.
The maximum error is 1.30051e-05, while requested accuracy is 1e-06.

The solution was obtained on a mesh of 1500 points. The maximum error is 1.301e-05. There were 71957 calls to the ODE function. There were 59 calls to the BC function. The parameter is 1.095 The final omega value is 1.066

Warning: Ignoring extra legend entries.

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

The solution was obtained on a mesh of 1500 points. The maximum error is 9.985e-06. There were 76180 calls to the ODE function. There were 60 calls to the BC function. The parameter is 1.184 The final omega value is 1.179

%Program to calculate the inital data

%Considering a spherically symmetric scalar field

%\psi(t,r)=\psi_0(r)exp(i\omega t)

%Line element in polar-areal gauge

%ds^2=-\alpha(r)^2 dt^2+A(r)^2 dr^2+r^2d\Omega^2

function onestar2

clear;

clc;

format long;

infinity = 15;

omega=0.50;

x_init=1e-5:0.01:infinity;

global phi_c

for phi_c=[0.07 0.04 0.02 0.01]

%solinit = bvpinit(linspace(1e-5,infinity,1000),[1 1 phi_c 0],omega);

solinit = bvpinit(x_init,[1 1 phi_c 0],omega);

options = bvpset('stats','on','RelTol',1e-6);

%condition=true;

%while condition

sol = bvp5c(@bsode,@bsbc,solinit,options);

%condition = sol.stats.maxerr >= 1e-5;

%end

r_data = sol.x;

f = sol.y;

i=1;

while (f(3,i)>1e-5)

i=i+1;

end

r_data = r_data(1:i);

%Calculating ADM Mass

A1=f(1,1:i);

mass=(r_data./2).*(1-A1.^(-2));

%Scaling alpha and omega

al_last = 1./f(2,i);

al_calc = al_last.*f(2,1:i);

omega_f = al_last*sol.parameters;

fprintf('\n');

fprintf('The parameter is %7.3f \n',...

sol.parameters)

fprintf('The final omega value is %7.3f \n',...

omega_f)

hold on;

figure(1)

plot(r_data,mass);

axis([0 infinity 0 0.7]);

title('Mass')

xlabel('r')

ylabel('M')

legend('0.07', '0.04', '0.02', '0.01')

hold on;

figure(2)

plot(r_data,al_calc);

axis([0 infinity 0.7 1.1]);

title('Lapse')

xlabel('r')

ylabel('\alpha')

legend('0.07', '0.04', '0.02', '0.01')

hold on;

figure(3)

plot(r_data,f(3,1:i));

axis([0 infinity 0 0.1]);

title('ScalarField')

xlabel('r')

ylabel('\Psi_0')

legend('0.07', '0.04', '0.02', '0.01')

hold off

% figure(4)

% plot(r_data,f(4,:));

% axis([0 infinity 0 2]);

% title('FieldDerivative')

% xlabel('r')

% ylabel('Q')

end

% --------------------------------------------------------------------------

%Function to decalre the complete set of differential equations

%y=[y(1),y(2),y(3),y(4)]=[A,alpha,psi_0,Q]

%V=1/2 m^2 \psi\psi* (free field term)

%dV/d\psi=1/2 m^2 \psi*

function dfdr = bsode(r,y,p)

m = 1;

%omega = 0.97;

dfdr = [ (1/2)*((y(1)/r)*(1-y(1)^2)+4*pi*r*y(1)*(y(3)^2*y(1)^2*(m^2+p(1)^2/y(2)^2)+y(4)^2))

(y(2)/2)*((y(1)^2-1)/r+4*pi*r*(y(3)^2*y(1)^2*(p(1)^2/y(2)^2-m^2)+y(4)^2))

y(4)

-(1+y(1)^2-4*pi*r^2*y(3)^2*y(1)^2*m^2)*(y(4)/r)-(p(1)^2/y(2)^2-m^2)*y(3)*y(1)^2];

% --------------------------------------------------------------------------

%Function file for declaring the boundary conditions

%Regularity at origin : A(0)=1, \psi_0(0)=\psi_c, Q(0)=0

%Asymptotic flat condition: A(r->\inf)=1, \alpha(r->inf)=1, \psi_0(r->inf)=0

end

function res = bsbc(ya,yb,p)

global phi_c

res = [ya(1)-1 %fixed

%yb(1)-1

ya(2)-1 %fixed

%yb(2)-1

ya(3)-phi_c %fixed

yb(3) %fixed

ya(4) ]; %fixed

end

6 个评论
显示 4更早的评论隐藏 4更早的评论

T 2024-9-3

编辑：Torsten 2024-9-4

在 MATLAB Online 中打开

I have scaled the equations and wrote this code keeping in mind the above code. Can you tell what am i doing wrong?

%Program to calculate the inital data
%Considering a spherically symmetric scalar field
%\psi(t,r)=\psi_0(r)exp(i\omega t)
%Line element in polar-areal gauge
%ds^2=-\alpha(r)^2 dt^2+A(r)^2 dr^2+r^2d\Omega^2
onestar2()
Error using bvp4c (line 196)
Unable to solve the collocation equations -- a singular Jacobian encountered.

Error in solution>onestar2 (line 30)
    sol = bvp4c(@bsode,@bsbc,solinit,options);
function onestar2
clear;
clc;
format long;
infinity = 15;
omega=0.50;
x_init=1e-5:0.01:infinity;
global phi_c
for phi_c=[0.07 0.04 0.02 0.01]
%solinit = bvpinit(linspace(1e-5,infinity,1000),[1 1 phi_c 0],omega);
    solinit = bvpinit(x_init,[1 1 phi_c 0],omega);
    options = bvpset('stats','on','RelTol',1e-6);
    %condition=true;
    %while condition
    
    sol = bvp4c(@bsode,@bsbc,solinit,options);
        %condition = sol.stats.maxerr >= 1e-5;
    %end
    r_data = sol.x;
    f = sol.y;
    
    i=1;
    
    while (f(3,i)>1e-5)
        i=i+1;
    end
    
    r_data = r_data(1:i);
    %Calculating ADM Mass
    A1=f(1,1:i);
    mass=(r_data./2).*(1-A1.^(-2));
    %Scaling alpha and omega
    al_last = 1./f(2,i);
    al_calc = al_last.*f(2,1:i); 
    omega_f = al_last*sol.parameters;
    fprintf('\n');
    fprintf('The parameter is %7.3f \n',...
        sol.parameters)
    fprintf('The final omega value is %7.3f \n',...
        omega_f)
    
    hold on;
    figure(1)
    plot(r_data,mass);
    axis([0 infinity 0 0.7]);
    title('Mass')
    xlabel('r')
    ylabel('M')
    legend('0.07', '0.04', '0.02', '0.01')
    hold on;
    
    figure(2)
    plot(r_data,al_calc);
    axis([0 infinity 0.7 1.1]);
    title('Lapse')
    xlabel('r')
    ylabel('\alpha')
    legend('0.07', '0.04', '0.02', '0.01')
    hold on;
    
    figure(3)
    plot(r_data,f(3,1:i));
    axis([0 infinity 0 0.1]);
    title('ScalarField')
    xlabel('r')
    ylabel('\Psi_0')
    legend('0.07', '0.04', '0.02', '0.01')
hold off
% figure(4)
% plot(r_data,f(4,:));
% axis([0 infinity 0 2]);
% title('FieldDerivative')
% xlabel('r')
% ylabel('Q')
end
end
% --------------------------------------------------------------------------
%Function to decalre the complete set of differential equations
%y=[y(1),y(2),y(3),y(4)]=[A,alpha,psi_0,Q]
%V=1/2 m^2 \psi\psi* (free field term)
%dV/d\psi=1/2 m^2 \psi*
function dfdr = bsode(r1,y,p)
   
m = 1;
G = 1;
%omega = 0.97;
a=y(1);
alpha1 = y(2);
psi0_1 = y(3);
phi = y(4);
dadr1 = (1/2) * a * ((1 - a^2)/r1 + 4 * pi * G * r1 * (psi0_1^2 * a^2 * (1 + 1 / alpha1^2)/(4*pi*G) + phi^2));
dalphadr1 = (alpha1/2) * (((a^2 - 1)/r1 + 4 * pi * G * r1 * (psi0_1^2 * a^2 * (1 / alpha1^2 - 1)/(4*pi*G) + phi^2)));
dpsidr1 = phi;
dphidr1 = -(1 + a^2 - r1^2 * psi0_1^2 * a^2) * (phi/r1) - (1 / alpha1^2 - 1)/(4*pi*G)^0.5 * psi0_1 * a^2;
dfdr = [dadr1;dalphadr1;dpsidr1;dphidr1];
%a , alpha1, psi0_1, phi = y
%dadr1 = (1/2) * a * ((1 - a^2)/r1 + 4 * np.pi * G * r1 * (psi0_1^2 * a^2 * (1 + 1 / alpha1^2)/(4*np.pi*G) + phi^2))
%dalphadr1 = (alpha1/2) * (((a^2 - 1)/r1 + 4 * np.pi * G * r1 * (psi0_1^2 * a^2 * (1 / alpha1^2 - 1)/(4*np.pi*G) + phi^2)))
%dpsidr1 = phi
%dphidr1 = -(1 + a^2 - r1^2 * psi0_1^2 * a^2) * (phi/r1) - (1 / alpha1^2 - 1)/(4*np.pi*G)^0.5 * psi0_1 * a^2
end
% --------------------------------------------------------------------------
%Function file for declaring the boundary conditions
%Regularity at origin : A(0)=1, \psi_0(0)=\psi_c, Q(0)=0
%Asymptotic flat condition: A(r->\inf)=1, \alpha(r->inf)=1, \psi_0(r->inf)=0 
function res = bsbc(ya,yb,p)
global phi_c
res = [ya(1)-1 %fixed
       %yb(1)-1 
       ya(2)-1 %fixed
       %yb(2)-1 
       ya(3)-phi_c %fixed
       yb(3)    %fixed
       ya(4) ]; %fixed
end

Torsten 2024-9-4

编辑：Torsten 2024-9-4

I removed the syntax errors, but MATLAB is not able to solve your equations (see above).

Since you fixed the values of all variables at x = 0, you should try an initial-value solver like ode45 to see what you get for different values of p. But I doubt that fixing all variables at x = 0 is the correct choice of the boundary conditions because it seems you have a singular point there.

T 2024-9-18

编辑：T 2024-9-18

the code i provided in the question is for mass-radius plot for a single boson star for different phi_c values. i want to plot the mass-radius plot for a system of boson stars for different range of phi_c values where each point of the plot denotes one star with one phi_c value. can i do that? if yes, how. i want a single curve which represents the mass-radius plot for system of stars. the final should resemble the blue line. the code should work in such a way that the odes are being solved with respect to the boundary conditions over the range of the phi_c given and finally it plots the final maximum mass and corresponding radius obtained for each phi_c value after which the graph gets stable for one star and plot it in one curve.

请先登录，再进行评论。

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

Solving bvp using shooting method to solve coupled first order odes

6 个评论
显示 4更早的评论隐藏 4更早的评论

回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

Solving bvp using shooting method to solve coupled first order odes

6 个评论 显示 4更早的评论隐藏 4更早的评论

回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

6 个评论
显示 4更早的评论隐藏 4更早的评论