如何在Scilab中求解带导数边界条件的含y'平方项的非线性二阶微分方程

待求解模型 #

给定常微分方程与边界条件：

-cos(y)y''+sin(y)y'^2+sin(y)=0, y=y(t), 边界条件为 y'(0)=0, y'(1)=0

该问题为带y'^2非线性项的两点边值问题，可通过Scilab内置边值求解器完成数值计算。

求解前置推导 #

Scilab的边值求解器仅支持一阶方程组形式的模型，因此需要先对原二阶方程做降阶处理：

• 定义状态变量：令y1(t) = y(t)（原未知函数），y2(t) = y'(t)（原函数的一阶导数）
• 状态导数推导：
  1. 第一个状态的导数：y1'(t) = y2(t)
  2. 从原方程移项解出二阶导y''(t)（即y2'(t)）：
     代入状态变量后原方程变为-cos(y1)y2' + sin(y1)y2^2 + sin(y1) = 0，整理得：
     y2'(t) = (sin(y1) * y2^2 + sin(y1)) / cos(y1) = tan(y1)*(y2^2 + 1)
• 边界条件转写：原边界条件对应状态变量的约束为t=0时y2(0)=0，t=1时y2(1)=0，共2个边界约束，匹配二阶方程的求解自由度。

代码实现步骤 #

• 第一步：定义一阶方程组的右端计算函数
  加极小量做分母保护，避免cos(y1)接近0时出现数值奇异：

  function dydt = odefun(t, yvec)
      y1 = yvec(1);
      y2 = yvec(2);
      eps_val = 1e-8;
      dydt(1) = y2;
      dydt(2) = (sin(y1) * y2^2 + sin(y1)) / (cos(y1) + eps_val);
  endfunction
  

• 第二步：定义边界残差函数
  求解器要求左右两端点的残差为0，对应设置两个边界约束：

  function res = bcfun(yL, yR)
      // yL为t=0处的状态值，yR为t=1处的状态值
      res(1) = yL(2); // 左边界：y'(0) = 0
      res(2) = yR(2); // 右边界：y'(1) = 0
  endfunction
  

• 第三步：构造初始猜测解
  非线性边值问题的求解收敛高度依赖初猜合理性。注意该方程存在平凡常数解：当y(t) = n*pi（n为整数）时，y'=0、y''=0，代入原方程可验证完全满足等式与边界条件。因此初猜可设置为接近这类常数解的形式，避免求解发散：

  // 在[0,1]区间生成20个均匀节点作为初始计算网格
  t_mesh = linspace(0, 1, 20);
  // 初猜y1取pi/4（靠近0常数解），y2全为0
  y_init = [pi/4 * ones(1, length(t_mesh)); zeros(1, length(t_mesh))];
  sol_init = bvpinit(t_mesh, y_init);
  

• 第四步：调用边值求解器计算

  sol = bvp4c(odefun, bcfun, sol_init);
  

• 第五步：结果验证与可视化
  对求解结果做插值，绘制原函数与一阶导数曲线，直观验证边界条件是否满足：

  t_dense = linspace(0, 1, 100);
  y_dense = bvpval(sol, t_dense);
  plot(t_dense, y_dense(1,:), 'r', t_dense, y_dense(2,:), 'b--');
  legend('y(t)', "y'(t)");
  xlabel('t');
  

常见问题排查 #

• 求解不收敛、报错发散的处理：
  • 拉近初猜值与已知平凡解的距离，不要选择离n*pi太远的初猜起点
  • 适当加密初始计算网格，降低初始插值带来的求解误差
  • 检查代码中分母保护的epsilon取值，避免步长过大时跳过奇异点
• 结果不满足边界条件的处理：
  • 核对边界残差函数的左右端点逻辑，避免将yL和yR的约束写反
  • 调小求解器的相对误差容限，提高数值计算精度

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