利用MATLAB数值求解边值问题——Berman问题

在偏微分方程的课堂上学习了Berman边界值问题

为了研究U235的扩散、气体的冷凝，Berman提出以下模型：

利用摄动方法，将非线性问题转化为线性问题，得出如下数学模型：

利用bvp4c函数数值求解，代码如下：

function bvp_example
    R_values = [25.445, 14.678, 13.226]; %不同的雷诺数
    colors = ['r', 'g', 'b']; 
    solinit = bvpinit(linspace(0, 1, 10), [-2 -90 -3 -10]); %赋给的初始值
    figure;
    hold on;
    xlabel('x');
    ylabel('f''(x)');
    title('f''图像 (不同R值)');
    grid on;
    for i = 1:length(R_values) %对不同雷诺数分别画图
        R = R_values(i);
        sol = bvp4c(@(x, y) odefun(x, y, R), @bcfun, solinit); %bvp4c函数
        x = linspace(0, 1, 100);
        y = deval(sol, x);
        plot(x, y(2, :), colors(i), 'DisplayName', ['R = ', num2str(R)]);
    end
    legend('show');
    function dydx = odefun(x, y, R) %将高阶常微分方程转化为一阶常微分方程组
        f = y(1);
        fp = y(2);
        fpp = y(3);
        fppp = y(4);
        dydx = [fp; fpp; fppp; R * (f * fppp - fp * fpp)];
    end
    function res = bcfun(ya, yb) %边界值
        res = [ya(1); ya(3); yb(2); yb(1) - 1];
    end
end

代码运行的结果如下图：