带约束非线性边值问题求解：MATLAB数值实现优化问询

问题描述 #

我需要求解如下带约束的非线性边值问题（对应方程见图片），其中S、k_i、x_f、α₁、θ为已知参数，目标求解h(x)、p和θ_d。

我的求解思路：采用有限差分法构建数值求解方案——对常微分方程使用二阶中心差分得到差分方程，结合Dirichlet边界条件h₀=h_N=0，通过一阶中心差分处理接触角边界条件得到虚点方程，以此补充求解p和θ_d所需的方程。

现欲在MATLAB中实现迭代求解方案：以S=0时的解析解作为初始猜测，每次迭代需求解非线性方程组。我考虑使用fsolve，但手动为每个h_i定义函数句柄在网格点数较多时过于繁琐（仅N=4时可手动实现）。请问：

1. 是否有更简便的实现方法？
2. 我的求解思路是否合理？

附N=4时的MATLAB测试代码：

clear; clc;
% Set problem parameters
xf = 6.17415;
theta = 39*(pi/180);
S = 1; ki=.37; alpha1 = 2.69;
% Set numerical parameters
L = 2*xf;
n = 4;
dx = L/n;
x=-xf:dx:xf;
iters = 5; % number of iterations

uold = zeros(1,n+3);   % First n-1 elements are h_1, ..., h_{n-1}; Last 2 elements are p and \theta_d
hold = (1 + (sinh(x-xf)-sinh(x+xf))/sinh(2*xf))*.809827;
uold(1,1:n+1) = hold;
uold(1,n+2:n+3) = [.809827,0];
plot(x,hold)
hold on
F = @(y) [1/(dx^2+3)*y(2)-(dx^2+2)/(dx^2+3)*y(1)-S/(dx^2+3)*exp(-2*ki*((x(2)+xf)+alpha1*hold(1,2)))+y(4)/(dx^2+3);
    1/(dx^2+3)*(y(3)+y(1))-(dx^2+2)/(dx^2+3)*y(2)-S/(dx^2+3)*exp(-2*ki*((x(3)+xf)+alpha1*hold(1,3)))+y(4)/(dx^2+3);
    1/(dx^2+3)*y(2)-(dx^2+2)/(dx^2+3)*y(3)-S/(dx^2+3)*exp(-2*ki*((x(4)+xf)+alpha1*hold(1,4)))+y(4)/(dx^2+3);
    y(4)+2*y(1)-dx*tan(theta-y(5))-S;
    y(5)+theta-atan(2/dx*y(3)-S/dx*exp(-4*ki*xf)+y(4)/dx);];

for j=1:iters
    u = fsolve(F,uold);
    plot(x,u(1,1:n+1))
    disp(u(1,n+2:n+3))
    uold = u;
end
hold off

解答 #

一、求解思路合理性 #

你的求解思路完全合理：

• 二阶中心差分对光滑解的精度足够，是这类非线性常微分方程离散的标准选择；
• 用虚点法处理接触角边界条件（将Neumann类边界转化为代数方程）是边值问题的常用技巧，能有效补充自由度，与Dirichlet边界条件配合构建完整的方程组；
• 以S=0时的解析解作为初始猜测，能大幅降低fsolve的收敛难度，避免迭代陷入局部极值。

二、简化fsolve函数句柄的实现方法 #

手动枚举每个差分方程的方式确实不适用于大网格，推荐通过循环+函数封装动态构建残差向量，无需手动定义每个h_i的方程：

1. 核心优化思路

将未知量向量y拆分为三部分：

• y(1:n+1)：所有网格点的h值（包含边界点h₀=h_N=0）；
• y(n+2)：参数p；
• y(n+3)：参数θ_d。

通过循环遍历内部网格点，自动生成每个点的差分残差，再添加接触角边界条件对应的残差，最终组合成完整的残差向量。

2. 优化后的代码示例

clear; clc;
% 问题参数
xf = 6.17415;
theta = 39*(pi/180);
S = 1; ki=.37; alpha1 = 2.69;
% 数值参数（可任意调整网格点数）
n = 20;
L = 2*xf;
dx = L/n;
x = -xf:dx:xf;
iters = 5;

% 初始猜测：S=0时的解析解
h_init = (1 + (sinh(x-xf)-sinh(x+xf))/sinh(2*xf))*.809827;
uold = [h_init, .809827, 0]; % 前n+1个是h，最后两个是p、theta_d
plot(x, h_init)
hold on

% 定义残差函数，动态生成方程组
F = @(y) compute_residual(y, x, dx, S, ki, alpha1, theta, n);

for j=1:iters
    u = fsolve(F, uold);
    plot(x, u(1:n+1))
    disp(['迭代', num2str(j), '：p=', num2str(u(n+2)), '，theta_d=', num2str(u(n+3))])
    uold = u;
end
hold off

% 残差计算子函数
function res = compute_residual(y, x, dx, S, ki, alpha1, theta, n)
    res = zeros(n+2, 1); % n-1个内部点差分方程 + 2个边界条件方程
    h = y(1:n+1);
    p = y(n+2);
    theta_d = y(n+3);
    C = 1/(dx^2 + 3); % 统一系数因子
    
    % 生成内部点的差分残差（i=2到n对应x(2)到x(n)，即内部网格点）
    for i=2:n
        if i == 2
            % 左内部点，邻点为i=1（边界h0=0）和i=3
            res(i-1) = C*h(i+1) - C*(dx^2 + 2)*h(i) - C*S*exp(-2*ki*((x(i)+xf)+alpha1*h(i))) + C*p;
        elseif i == n
            % 右内部点，邻点为i=n-1和i=n+1（边界hN=0）
            res(i-1) = C*h(i-1) - C*(dx^2 + 2)*h(i) - C*S*exp(-2*ki*((x(i)+xf)+alpha1*h(i))) + C*p;
        else
            % 中间内部点，二阶中心差分
            res(i-1) = C*(h(i+1)+h(i-1)) - C*(dx^2 + 2)*h(i) - C*S*exp(-2*ki*((x(i)+xf)+alpha1*h(i))) + C*p;
        end
    end
    
    % 添加接触角边界条件对应的残差
    res(n) = p + 2*h(2) - dx*tan(theta - theta_d) - S;
    res(n+1) = theta_d + theta - atan(2/dx*h(n) - S/dx*exp(-4*ki*xf) + p/dx);
end

3. 关键优化点

• 将残差计算封装为独立函数，通过循环自动生成所有内部点的差分方程，无论网格点数n多大都无需手动修改；
• 未知量向量结构统一，fsolve可直接调用，代码可读性和可维护性大幅提升；
• 后续调整差分格式或边界条件，仅需修改循环内的逻辑即可。

额外建议 #

• 若网格点数极大，可考虑用数组切片替代循环进一步提升计算效率（如用h(3:n+1)+h(1:n-1)处理中间点的邻点和）；
• 可通过optimset('Display','iter')给fsolve设置迭代显示选项，方便调试收敛情况；
• 若S偏离0较大，可采用延续法：从S=0开始逐步增大S进行求解，用前一步的解作为后一步的初始猜测，进一步保证收敛。

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