风落平川 2024-05-22 12:37 采纳率: 96.8%
浏览 19
已结题

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

请用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为剖分节点上的数值解

img

  • 写回答

20条回答 默认 最新

  • 专家-赤兔[在线] 全栈领域优质创作者 2024-05-22 12:37
    关注

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

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

    1. 定义系数函数 p(x)q(x)r(x)
    2. 计算边界条件的值 alphabeta
    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
    

    接下来,计算边界条件的值 alphabeta。题目已经给出 y(a) = alpha = 0y(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 值,解会更接近于实际的连续解。

    展开全部

    本回答被题主选为最佳回答 , 对您是否有帮助呢?
    评论 编辑记录
    专家-赤兔[在线] 回复 风落平川 2024-05-22 12:38

    要使用有限差分法求解这个一维边值问题,我们首先需要定义三个辅助函数:p_func(x), q_func(x), 和 r_func(x) 分别对应方程中的系数 p(x), q(x), 以及右端函数 r(x)。然后,我们将根据这些函数和给定的边界条件在区间 [a, b] 上使用中心差分方法近似二阶导数,并构造一个线性系统来求解数值解。

    以下是具体的步骤和代码实现:

    辅助函数定义

    首先,我们定义方程系数和右端函数的表达式。

    function pd = p_func(x)
    pd = 3; % 因为 p(x) = 3 是常数
    end
    
    functionqd = q_func(x)
    qd = 2; % 因为 q(x) = 2 是常数
    end
    
    function rd = r_func(x)
    rd = sin(x) + 3*cos(x); % 右端函数 r(x)
    end
    

    主函数实现

    接下来,我们实现主函数 myFD,使用上述辅助函数计算数值解。

    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] 剖分的等分数
    
    % 参数初始化
    h = (b - a) / N; % 空间步长
    x = a:h:b; % 剖分节点
    y = zeros(1, N+1); % 初始化y的数值解数组
    
    % 差分矩阵构建
    A = spdiags([-1 2 -1], -1:1, N+1, N+1); % 中心差分二阶导数矩阵
    I = eye(N+1); % 单位矩阵
    
    % 方程系数矩阵
    coeff_p = -p_func(x(2:end-1)) ./ (2*h);
    coeff_q = q_func(x(2:end-1));
    
    % 构造线性系统
    A_coeff = I + coeff_p * diag(ones(1,N-1),-1) + coeff_p * diag(ones(1,N-1),1) + coeff_q * I;
    rhs = r_func(x(2:end-1)) * h^2;
    
    % 添加边界条件
    A_coeff(1,:) = [1, zeros(1,N)]; % y(a)=alpha
    A_coeff(end,:) = [zeros(1,N), 1]; % y(b)=beta
    rhs([1,end]) = [alpha, beta];
    
    % 解线性系统
    y(2:end-1) = A_coeff \ rhs; % 解得内部节点的y值
    
    end
    

    使用示例

    现在,我们可以调用 myFD 函数求解算例1:

    alpha = 0; % 左边界条件
    beta = 1; % 右边界条件
    a = 0; % 区间左端点
    b = pi/2; % 区间右端点
    N = 50; % 剖分次数
    
    [x, y] = myFD(alpha, beta, a, b, N);
    

    这将返回区间 [a, b] 上的剖分节点 x 和对应的数值解 y。请注意,由于有限差分法的近似性质,数值解与精确解 y(x) = sin(x) 在边界和内部节点上可能会有微小误差。为了提高精度,可以增加剖分次数 N

    回复
查看更多回答(19条)
编辑
预览

报告相同问题?

问题事件

  • 系统已结题 6月5日
  • 已采纳回答 5月29日
  • 创建了问题 5月22日

悬赏问题

  • ¥15 PADS Logic 原理图
  • ¥15 PADS Logic 图标
  • ¥15 电脑和power bi环境都是英文如何将日期层次结构转换成英文
  • ¥20 气象站点数据求取中~
  • ¥15 如何获取APP内弹出的网址链接
  • ¥15 wifi 图标不见了 不知道怎么办 上不了网 变成小地球了
手机看
程序员都在用的中文IT技术交流社区

程序员都在用的中文IT技术交流社区

专业的中文 IT 技术社区,与千万技术人共成长

专业的中文 IT 技术社区,与千万技术人共成长

关注【CSDN】视频号,行业资讯、技术分享精彩不断,直播好礼送不停!

关注【CSDN】视频号,行业资讯、技术分享精彩不断,直播好礼送不停!

客服 返回
顶部