This example shows how to solve Emden's equation, which is a boundary value problem with a singular term that arises in modeling a spherical body of gas.
After reducing the PDE of the model using symmetry, the equation becomes a second-order ODE defined on the interval ,
At , the term is singular, but symmetry implies the boundary condition . With this boundary condition, the term is well defined as . For the boundary condition , the BVP has an analytical solution that you can compare to a numeric solution calculated in MATLAB®,
The BVP solver
bvp4c can solve singular BVPs that have the form
The matrix must be constant and the boundary conditions at must be consistent with the necessary condition . Use the
'SingularTerm' option of
bvpset to pass the
S matrix to the solver.
You can rewrite Emden's equation as a system of first-order equations using and as
The boundary conditions become
In the required matrix form, the system is
In matrix form it is clear that and .
To solve this system of equations in MATLAB, you need to code the equations, boundary conditions, and options before calling the boundary value problem solver
bvp4c. You either can include the required functions as local functions at the end of a file (as done here), or save them as separate, named files in a directory on the MATLAB path.
Create a function to code the equations for . This function should have the signature
dydx = emdenode(x,y), where:
x is the independent variable.
y is the dependent variable.
dydx(1) gives the equation for , and
dydx(2) the equation for .
These inputs are automatically passed to the function by the solver, but the variable names determine how you code the equations. In this case:
function dydx = emdenode(x,y) dydx = [y(2) -y(1)^5]; end
The term that contains
S is passed to the solver using options, so that term is not included in the function.
Now, write a function that returns the residual value of the boundary conditions at the boundary points. This function should have the signature
res = emdenbc(ya,yb), where:
ya is the value of the boundary condition at the beginning of the interval.
yb is the value of the boundary condition at the end of the interval.
For this problem, one of the boundary conditions is for , and the other is for . To calculate the residual values, you need to put the boundary conditions into the form .
In this form the boundary conditions are
The corresponding function is then
function res = emdenbc(ya,yb) res = [ya(2) yb(1) - sqrt(3)/2]; end
Lastly, create an initial guess of the solution. For this problem, use a constant initial guess that satisfies the boundary conditions, and a simple mesh of five points between 0 and 1. Using many mesh points is unnecessary since the BVP solver adapts these points during the solution process.
guess = [sqrt(3)/2; 0]; xmesh = linspace(0,1,5); solinit = bvpinit(xmesh, guess);
Create a matrix for
S and pass it to
bvpset as the value of the
'SingularTerm' option. Finally, call
bvp4c with the ODE function, boundary condition function, initial guess, and option structure.
S = [0 0; 0 -2]; options = bvpset('SingularTerm',S); sol = bvp4c(@emdenode, @emdenbc, solinit, options);
Calculate the value of the analytical solution in .
x = linspace(0,1); truy = 1 ./ sqrt(1 + (x.^2)/3);
Plot the analytical solution and the solution calculated by
bvp4c for comparison.
plot(x,truy,sol.x,sol.y(1,:),'ro'); title('Emden Problem -- BVP with Singular Term.') legend('Analytical','Computed'); xlabel('x'); ylabel('solution y');
Listed here are the local helper functions that the BVP solver
bvp4c calls to calculate the solution. Alternatively, you can save these functions as their own files in a directory on the MATLAB path.
function dydx = emdenode(x,y) % equation being solved dydx = [y(2) -y(1)^5]; end %------------------------------------------- function res = emdenbc(ya,yb) % boundary conditions res = [ya(2) yb(1) - sqrt(3)/2]; end %-------------------------------------------