Matlab求解边值问题：如何将符号解转为纯数值解？

我正在编写Matlab脚本，用指定公式求解边值问题（BVP）。当前代码能运行，但输出的解向量包含y2、y3等符号变量，希望得到全数值的解向量，该如何实现？

原代码 #

%Valors de la frontera:  y''(ti)+p(ti)·y(ti)+q(ti)·y(ti)=f(ti) %  
p=@(t,y) -t;
q=@(t,y) 2;
f=@(t,y) 0;
n=4;
a=0;
b=1;
h=(b-a)/n;
y0=0;
yf=1;
Y=[y0];
% Creamos el vector Y dinámicamente
K = sym('y', [1, n-1]);
Y=[Y K yf];
T = a + (0:n) * h;
t=a;
y=y0;

for i = 2:n
    t = T(i);
    y_prev = Y(i-1);
    y_current = Y(i);
    y_next = Y(i+1);
    equation = (1 - (h/2) * p(t, y_current)) * y_prev + (-2 + h^2 * q(t, y_current)) * y_current + (1 + (h/2) * p(t, y_current)) * y_next - h^2 * f(t, y_current);
    Y(i) = solve(equation, y_current);
end

Y

当前输出 #

Y = [0, (31*y2)/60, (900*y3)/1273, 36917/44880, 1]

问题原因与解决方法 #

原代码错误地用符号变量逐个求解单个方程，导致变量互相依赖无法得到数值解。实际上，这是一个线性方程组求解问题，应该通过构建系数矩阵，用数值方法一次性求解所有内点的未知值。

具体步骤：

• 设中间n-1个点的函数值为未知向量x
• 根据差分格式，为每个内点建立线性方程，组成系数矩阵A和右端向量b
• 代入边界条件（y0=0，yf=1），求解线性方程组A*x = b
• 拼接边界值和求解结果，得到全数值的解向量

修改后的代码 #

% 边值问题：y''(t) + p(t)y'(t) + q(t)y(t) = f(t)，边界条件y(a)=y0, y(b)=yf
p=@(t) -t;  % 原代码中p的y参数未使用，简化为仅关于t的函数
q=@(t) 2;
f=@(t) 0;
n=4;
a=0;
b=1;
h=(b-a)/n;
y0=0;
yf=1;

% 生成节点
T = a + (0:n)*h;
% 内点数量：n-1个
m = n-1;
% 初始化系数矩阵A和右端向量b
A = zeros(m,m);
b = zeros(m,1);

% 为每个内点构建方程
for i = 1:m
    t = T(i+1);  % 第i个内点对应的t值
    % 差分格式的系数
    coeff_prev = 1 - (h/2)*p(t);
    coeff_curr = -2 + h^2*q(t);
    coeff_next = 1 + (h/2)*p(t);
    
    % 处理第一个内点（前一个点是边界y0）
    if i == 1
        A(i,i) = coeff_curr;
        A(i,i+1) = coeff_next;
        b(i) = h^2*f(t) - coeff_prev*y0;
    % 处理最后一个内点（后一个点是边界yf）
    elseif i == m
        A(i,i-1) = coeff_prev;
        A(i,i) = coeff_curr;
        b(i) = h^2*f(t) - coeff_next*yf;
    % 中间内点
    else
        A(i,i-1) = coeff_prev;
        A(i,i) = coeff_curr;
        A(i,i+1) = coeff_next;
        b(i) = h^2*f(t);
    end
end

% 求解线性方程组
x = A\b;
% 拼接边界值和内点解，得到完整的解向量Y
Y = [y0; x; yf];
% 转换为行向量输出（可选）
Y = Y';

disp('全数值解向量Y：');
disp(Y);

输出结果 #

全数值解向量Y：
         0    0.1978    0.4898    0.8226    1.0000

内容的提问来源于stack exchange，提问作者Marc Vilà Llinàs