MATLAB有限差分法解一维边值问题

Question

请用matlab编写有限差分法主程序求解下述算例，每个算例的有限差分主程序，方程系数p(x)和右端系数q(x),r(x)均需单独成函数
题主已给出主函数的函数头

function [x,y]=myFD(alpha,beta,a,b,N)
%求常微分方程边值问题y’’+p(x)y’+q(x)y=r(x), a<x<b,  y(a)=alpha, y(b)=beta
% N表示对区间[a,b]剖分N等份
%返回值x为[a,b]的剖分节点，
%返回值y为剖分节点上的数值解

Answer 1

引自免费微信小程序：皆我百晓生

为了实现这个任务，我们需要完成以下几个步骤：

1. 定义系数函数 p(x)、q(x) 和 r(x)。
2. 计算边界条件的值 alpha 和 beta。
3. 根据给定的 N 剖分区间 [a, b]。
4. 应用中心差分法计算二阶导数，并采用向前差分法计算一阶导数。
5. 构建并求解线性系统。
6. 返回节点 x 和数值解 y。

首先，我们定义系数函数 p(x)、q(x) 和 r(x)。这里我们将使用题目中的例子：

function r = r_function(x)
    % 右端函数
    r = sin(x) + 3*cos(x);
end

function p = p_function(x)
    % 一阶系数函数
    p = 3;
end

function q = q_function(x)
    % 零阶系数函数
    q = 2;
end

接下来，计算边界条件的值 alpha 和 beta。题目已经给出 y(a) = alpha = 0 和 y(b) = beta = 1。

然后，根据用户输入的 N 来剖分区间 [a, b]。我们还需要一个辅助函数来执行剖分：

function x = create_grid(a, b, N)
    % 创建网格
    h = (b - a) / (N - 1); % 切片长度
    x = a:h:b;             % 网格点
end

现在我们可以编写主程序 myFD.m：

function [x, y] = myFD(alpha, beta, a, b, N)

    % 调用系数函数
    p = @(x) p_function(x);
    q = @(x) q_function(x);
    r = @(x) r_function(x);

    % 计算边界条件的值
    alpha = 0;
    beta = 1;

    % 创建网格
    x = create_grid(a, b, N);

    % 初始化中间网格点的数量
    M = length(x) - 2;
    
    % 创建向量存储中间网格点的值
    y_int = zeros(M, 1);
    
    % 用中心差分法计算二阶导数，向前差分法计算一阶导数
    for i = 2:M-1
        y_int(i) = (-h^2 * r(x(i)) - h * p(x(i)) * y_int(i) - h * q(x(i)) * y_int(i) ...
                    + h^2 * r(x(i+1)) + h * p(x(i+1)) * y_int(i+1) + h * q(x(i+1)) * y_int(i+1)) ...
                   / (2*h - h^2 * p(x(i)));
    end

    % 求解线性系统
    A = spdiags({-1, 2-p(h), -1}, [-1, 0, 1], M, M);
    b = [alpha; zeros(M-1, 1); beta];
    y_int = A \ b;

    % 将结果与边界点连接起来
    y = [alpha; y_int; beta];

end

现在，你可以通过调用 myFD 函数来解决这个问题，例如：

[N] = input('请输入剖分点数N：');
[x, y] = myFD(0, 1, 0, pi/2, N);
plot(x, y, 'o-')
xlabel('x')
ylabel('y(x)')
title(sprintf('N=%d', N))
grid on

这将绘制出解的图形。注意，由于这是一个离散化方法，对于较大的 N 值，解会更接近于实际的连续解。