本文说明如何使用 bvp4c 求解具有未知参数的边界值问题。
马蒂厄方程在区间 [0,π] 上定义为
y′′+(λ−2qcos(2x))y=0y'^{'} +(λ−2q cos(2x))y = 0y′′+(λ−2q cos(2x))y=0。
当参数 q=5 时，边界条件为
y′(0)=0，
y′(π)=0。
但这最多只能将 y(x) 确定为一个数乘，因此需要第三个条件来指定特定解，
y(0)=1。
要在 MATLAB 中对此方程组求解，您需要先编写方程组、边界条件和初始估计值的代码，然后再调用边界值问题求解器 bvp4c。您可以将所需的函数作为局部函数包含在文件末尾，或者将它们作为单独的命名文件保存在 MATLAB 路径上的目录中。
编写方程代码
创建一个函数以编写方程代码。此函数应具有签名 dydx = mat4ode(x,y,lambda)，其中：
x 是自变量。
y 是因变量。
lambda 是表示特征值的未知参数。
您可以用代换法 y1=y和y2=y′y_1=y 和 y_2=y′y1=y和y2=y′
将马蒂厄方程写成一阶方程组，
y1′=y2y_1' = y_2y1′=y2，
y2′=−(λ−2qcos(2x))y1y_2' = −(λ−2q cos(2x))y_1y2′=−(λ−2q cos(2x))y1。
则对应的函数是
function dydx = mat4ode(x,y,lambda) % equation being solveddydx = [y(2)-(lambda - 2*q*cos(2*x))*y(1)];end
注意：所有函数都作为局部函数包含在本文的末尾。
编写边界条件代码
现在，编写一个函数，该函数返回在边界点处的边界条件的残差值。此函数应具有签名 res = mat4bc(ya,yb,lambda)，其中：
ya 是在区间 [a,b] 开始处的边界条件的值。
yb 是在区间 [a,b] 结束处的边界条件的值。
lambda 是表示特征值的未知参数。
此问题在区间 [0,π] 内有三个边界条件。要计算残差值，您需要将边界条件设置为 g(x,y)=0 形式。在此形式中，边界条件是
y′(0)=0，
y′(π)=0，
y(0)−1=0。
则对应的函数是
function res = mat4bc(ya,yb,lambda) % boundary conditionsres = [ya(2)yb(2)ya(1)-1];end
创建初始估计值
最后，创建解的初始估计值。您必须对两个解分量 y1=y(x)y_1=y(x)y1=y(x)
和y2=y′(x)y_2=y'(x)y2=y′(x)以及未知参数 λ 提供初始估计值。
对于此问题，余弦函数满足三个边界条件，因此有助于提供较好的初始估计值。使用返回 y1y_1y1 和 y2y_2y2 的估计值的函数，编写 y 的初始估计值的代码。
function yinit = mat4init(x) % initial guess functionyinit = [cos(4*x)-4*sin(4*x)];end
使用区间为 [0,π] 的 10 点网格、初始估计值函数以及 λ 的估计值 15 调用 bvpinit。
lambda = 15;solinit = bvpinit(linspace(0,pi,10),@mat4init,lambda);
求解方程
使用 ODE 函数、边界条件函数和初始估计值调用 bvp4c。
sol = bvp4c(@mat4ode, @mat4bc, solinit);
参数值
打印 bvp4c 求得的未知参数 λ 的值。此值是马蒂厄方程的第四个特征值 (q=5)。
fprintf('Fourth eigenvalue is approximately %7.3f.\n',...sol.parameters)Fourth eigenvalue is approximately 17.097.
对解进行绘图
使用 deval 计算 bvp4c 在区间 [0,π] 中的 100 个点处计算的解。
xint = linspace(0,pi);Sxint = deval(sol,xint);
对两个解分量进行绘图。绘图显示了与第四个特征值 λ4=17.097λ_4=17.097λ4=17.097
相关联的特征函数（及其导数）。
plot(xint,Sxint)axis([0 pi -4 4])title('Eigenfunction of Mathieu''s Equation.') xlabel('x')ylabel('y')legend('y','y''')
局部函数
此处列出了 BVP 求解器 bvp4c 为计算解而调用的局部辅助函数。您也可以将这些函数作为它们自己的文件保存在 MATLAB 路径上的目录中。
function dydx = mat4ode(x,y,lambda) % equation being solvedq = 5;dydx = [y(2)-(lambda - 2*q*cos(2*x))*y(1)];end%-------------------------------------------function res = mat4bc(ya,yb,lambda) % boundary conditionsres = [ya(2)yb(2)ya(1)-1];end%-------------------------------------------function yinit = mat4init(x) % initial guess functionyinit = [cos(4*x)-4*sin(4*x)];end%-------------------------------------------
