Solve BVP with Unknown Parameter

This example shows how to use bvp4c to solve a boundary value problem with an unknown parameter.

Mathieu's equation is defined on the interval $[0, π]$ by

${y^{'}}^{'} + (λ - 2 q \cos (2 x)) y = 0$ .

When the parameter $q = 5$ , the boundary conditions are

$y^{'} (0) = 0$ ,

$y^{'} (π) = 0$ .

However, this only determines $y (x)$ up to a constant multiple, so a third condition is required to specify a particular solution,

$y (0) = 1$ .

To solve this system of equations in MATLAB®, you need to code the equations, boundary conditions, and initial guess before calling the boundary value problem solver bvp4c. You can either 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.

Code Equation

Create a function to code the equations. This function should have the signature dydx = mat4ode(x,y,lambda), where:

x is the independent variable.
y is the dependent variable.
lambda is an unknown parameter representing an eigenvalue.

You can write Mathieu's equation as a first-order system using the substitutions $y_{1} = y$ and $y_{2} = y^{'}$ ,

${y_{1}}^{'} = y_{2}$ ,

${y_{2}}^{'} = - (λ - 2 q \cos (2 x)) y_{1}$ .

The corresponding function is then

function dydx = mat4ode(x,y,lambda) % equation being solved
dydx = [y(2)
      -(lambda - 2*q*cos(2*x))*y(1)];
end

Note: All functions are included as local functions at the end of the example.

Code Boundary Conditions

Now, write a function that returns the residual value of the boundary conditions at the boundary points. This function should have the signature res = mat4bc(ya,yb,lambda), where:

ya is the value of the boundary condition at the beginning of the interval $[a, b]$ .
yb is the value of the boundary condition at the end of the interval $[a, b]$ .
lambda is an unknown parameter representing an eigenvalue.

This problem has three boundary conditions in the interval $[0, π]$ . To calculate the residual values, you need to put the boundary conditions into the form $g (x, y) = 0$ . In this form the boundary conditions are

$y^{'} (0) = 0$ ,

$y^{'} (π) = 0$ ,

$y (0) - 1 = 0$ .

The corresponding function is then

function res = mat4bc(ya,yb,lambda) % boundary conditions
res = [ya(2)
       yb(2)
       ya(1)-1];
end

Create Initial Guess

Lastly, create an initial guess of the solution. You must provide an initial guess for both solution components $y_{1} = y (x)$ and $y_{2} = y^{'} (x)$ , as well as the unknown parameter $λ$ . Only eigenvalues and eigenfunctions that are close to the initial guesses can be computed. To increase the likelihood that the computed eigenfunction corresponds to the fourth eigenvalue, you should choose an initial guess that has the correct qualitative behavior.

For this problem, a cosine function makes for a good initial guess since it satisfies the three boundary conditions. Code the initial guess for $y$ using a function that returns the guess for $y_{1}$ and $y_{2}$ .

function yinit = mat4init(x) % initial guess function
yinit = [cos(4*x)
        -4*sin(4*x)];
end

Call bvpinit using a mesh of 10 points in the interval $[0, π]$ , the initial guess function, and a guess of 15 for the value of $λ$ .

lambda = 15;
solinit = bvpinit(linspace(0,pi,10),@mat4init,lambda);

Solve Equation

Call bvp4c with the ODE function, boundary condition function, and initial guess.

sol = bvp4c(@mat4ode, @mat4bc, solinit);

Value of Parameter

Print the value of the unknown parameter $λ$ found by bvp4c. This value is the fourth eigenvalue ( $q = 5$ ) of Mathieu's equation.

fprintf('Fourth eigenvalue is approximately %7.3f.\n',...
            sol.parameters)

Fourth eigenvalue is approximately  17.097.

Plot Solution

Use deval to evaluate the solution computed by bvp4c at 100 points in the interval $[0, π]$ .

xint = linspace(0,pi);
Sxint = deval(sol,xint);

Plot both solution components. The plot shows the eigenfunction (and its derivative) associated with the fourth eigenvalue $λ_{4} = 17.097$ .

plot(xint,Sxint)
axis([0 pi -4 4])
title('Eigenfunction of Mathieu''s Equation.') 
xlabel('x')
ylabel('y')
legend('y','y''')

Figure contains an axes object. The axes object with title Eigenfunction of Mathieu's Equation., xlabel x, ylabel y contains 2 objects of type line. These objects represent y, y'.

Local Functions

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 = mat4ode(x,y,lambda) % equation being solved
q = 5;
dydx = [y(2)
      -(lambda - 2*q*cos(2*x))*y(1)];
end
%-------------------------------------------
function res = mat4bc(ya,yb,lambda) % boundary conditions
res = [ya(2)
       yb(2)
       ya(1)-1];
end
%-------------------------------------------
function yinit = mat4init(x) % initial guess function
yinit = [cos(4*x)
        -4*sin(4*x)];
end
%-------------------------------------------