Matlab bvp5c求解五阶微分方程边值问题的代码实现求助

问题背景 #

• 待求解的五阶微分方程形式如下：
  
• 求解首先需要明确边界条件，原始边界条件定义为：
  
• 由于需要补充额外条件实现方程系统闭合，首个边界条件需改写为其他形式，选取的改写形式为：
  
• 对新构建的问题做渐近求解，求解问题形式为：
  
• 该微分方程的解共对应5个特征根，其中2个为增长指数项，可通过如下条件消去：
  
• 基于上述条件推导得到2个边界条件：
  • 
  • 
• 此外补充两个边界条件：
  • 
  • F(0) = 1

现有实现 #

尝试使用Matlab内置的bvp5c工具实现该问题求解，但始终运行失败，已编写的代码如下：

eta0 = 20;
etamesh = [linspace(-eta0,0,100),linspace(0,eta0,100)];
Finit = [1; 0.5; 1; 0; 0];

solinit = bvpinit(etamesh,Finit)
sol = bvp5c(@fun_ode, @bc, solinit);

figure()
plot(sol.x,sol.y(1,:),'linewidth',1.5) % 绘制F曲线
ylim([-2 2]);
hold on
plot(sol.x,sol.y(3,:),'linewidth',1.5) % 绘制F''曲线
plot(sol.x,sol.y(5,:),'linewidth',1.5) % 绘制F''''曲线
grid on

legend("F","F''","F''''",'location','northwest')

function dFdeta = fun_ode(eta,F,region)
  dFdeta=[F(2);F(3);F(4);F(5);(1-F(1))/F(1)^3];
end

% 边界条件说明：
% F'''(-eta0) = F''''(-eta0) = 0
% F(0) = 1
% eta = eta0处[1] = [2] = 0

function res = bc(FL,FR)

alpha = cos(3*pi/5) + 1i*sin(3*pi/5);
beta  = cos(7*pi/5) + 1i*sin(7*pi/5);
gamma = abs(alpha);

  res = [FL(4,1) ; FL(5,1);...

         FR(1,1) - 1 ;...
         FR(1,1)-FL(1,2) ; FR(2,1)-FL(2,2) ; FR(3,1)-FL(3,2) ; FR(4,1)-FL(4,2) ;  FR(5,1)-FL(5,2) ;...
         
         (1-(beta+alpha))*FR(3,2) + (-(alpha+beta)+gamma^2)*FR(2,2) ; (-(beta+alpha)+gamma^2)*FR(3,2) + (FR(2,2))*gamma^2];
end   

现存问题 #

当前实现的核心难点有两点：

• 两个推导得到的耦合边界条件编写逻辑有误
• eta=0即F(0)处的解连续性条件设置不正确
  需要排查代码中的潜在错误，给出可运行的实现方案。

内容的提问来源于stack exchange，提问作者Wiss