Solving Singular BVPs :: Differential Equations (Mathematics)

Mathematics

Solving Singular BVPs

The function bvp4c solves a class of singular BVPs of the form

Equation 1: y prime = 1 / x times S times y + f(x,y).Equation 2: 0 = g(y(0),y(b)). (5-2)

It can also accommodate unknown parameters for problems of the form

Singular problems must be posed on an interval [0, b] with 0 --> b > . Use bvpset to pass the constant matrix to bvp4c as the value of the 'SingularTerm' integration property. Boundary conditions at x = 0 must be consistent with the necessary condition for a smooth solution, Sy(0) = 0 . An initial guess should also satisfy this necessary condition.

When you solve a singular BVP using

sol = bvp4c(@odefun,@bcfun,solinit,options)

bvp4c requires that your function odefun(x,y) return only the value of the f(x,y) term in Equation 5-2.

Example: Solving a BVP that Has a Singular Term

Emden's equation arises in modeling a spherical body of gas. The PDE of the model is reduced by symmetry to the ODE

on an interval [0, 1] . The coefficient 2 / x is singular at x = 0 , but symmetry implies the boundary condition y prime (0) = 0 . With this boundary condition, the term

is well-defined as approaches 0. For the boundary condition y(1) = square root of 3 divided by 2 , this BVP has the analytical solution

Note The demo emdenbvp contains the complete code for this example. The demo uses subfunctions to place all required functions in a single M-file. To run this example type emdenbvp at the command line. See BVP Solver Basic Syntax and Solving BVP Problems for more information.

Rewrite the problem as a first-order system and identify the singular term. Using a substitution and , write the differential equation as a system of two first-order equations

The boundary conditions become

Writing the ODE system in a vector-matrix form

the terms of Equation 5-2 are identified as

f(x,y) = a column vector containing two elements, first element is y sub 2, second element is minus y sub 1 to the fifth power

Code the ODE and boundary condition functions. Code the differential equation and the boundary conditions as functions that bvp4c can use.

function dydx = emdenode(x,y)
dydx = [  y(2) 
         -y(1)^5 ];
function res = emdenbc(ya,yb)
res = [ ya(2)
        yb(1) - sqrt(3)/2 ];

Setup integration properties. Use the matrix as the value of the 'SingularTerm' integration property.

S = [0,0;0,-2];
options = bvpset('SingularTerm',S);

Create an initial guess. This example starts with a mesh of five points and a constant guess for the solution.

Use bvpinit to form the guess structure

guess = [sqrt(3)/2;0];
solinit = bvpinit(linspace(0,1,5),guess);

Solve the problem. Use the standard bvp4c syntax to solve the problem.

sol = bvp4c(@emdenode,@emdenbc,solinit,options);

View the results. This problem has an analytical solution

The example evaluates the analytical solution at 100 equally-spaced points and plots it along with the numerical solution computed using bvp4c.

x = linspace(0,1);
truy = 1 ./ sqrt(1 + (x.^2)/3);
plot(x,truy,sol.x,sol.y(1,:),'ro');
title('Emden problem -- BVP with singular term.')
legend('Analytical','Computed');
xlabel('x');
ylabel('solution y');

Using Continuation to Make a Good Initial Guess Solving Multi-Point BVPs