数值计算讨论组

欢迎来到数值计算讨论组。 本论坛旨在增强师生交流的时效性,主要讨论数值计算相关问题,包括初等数值分析、矩阵计算、非线性方程(组)数值解法、微分方程数值解法等方面的内容。


    [微]中心差分格式求解两点边值问题 参考程序

    分享

    Admin
    Admin

    帖子数 : 46
    注册日期 : 10-05-13

    [微]中心差分格式求解两点边值问题 参考程序

    帖子 由 Admin 于 周三 十二月 08, 2010 9:06 pm

    本问题的难度比较低,主要是求解一个三对角线性方程组。该线性方程组的求解使用了函数TridiagAXb。下面的程序是主程序,后面会给出函数TridiagAXb的代码。
    代码:

    % 本程序使用中心差分格式计算两点边值问题
    % -u'' = f, u(0) = u(1) = 0
    % u_exact = 1 - (1-exp(-10))*x - exp(-10*x)
    % f = 100*exp(-10*x)
    %  n        误差结果-无穷范数 err_inf
    %  20        1.3742903775545e-2
    %  40        3.476685596693e-3
    %  80        8.712085690456e-4
    %  160      2.179740562154e-4
    % 2阶收敛速度
    %
    %  n        误差结果-无穷范数 err_L2
    %  20        9.277264457428e-3
    %  40        2.341306410922e-3
    %  80        5.867024245239e-4
    %  160      1.467616112220e-4
    % 2阶收敛速度


    clear; close all;format long;
    n = 160;
    h = 1/n;
    x = h:h:1-h;
    u_exact = zeros(n-1,1);
    u = zeros(n-1,1);
    b = zeros(n-1,1);
    for i = 1:n-1
        b(i) = f(x(i))*h^2;
        u_exact(i) = 1 - (1-exp(-10))*x(i) - exp(-10*x(i));
    end
    u = TridiagAXb(-1,2,-1,b);

    % 画图
    plot(x,u,'+',x,u_exact,'r')
    % 误差计算
    err_inf = max(abs(u - u_exact))
    err_L2 = 0;
    for i = 1:n-1
        err_L2 = err_L2 + h*(u(i) - u_exact(i))^2;
    end
    err_L2 = sqrt(err_L2)


    函数TridiagAXb

    代码:

    % 三对角方程求解函数  M*M
    % AX = B,  A的主对角线元素为b,右上次对角线元素为a,左下次对角线元素为c
    function X = TridiagAXb(a,b,c,B)
    M = length(B);
    X = zeros(M,1);
    D = zeros(1,M);
    D(1) = b;
    for k = 1:M-1
        temp = c/D(k);
        D(k+1) = b - a*temp;
        B(k+1) = B(k+1) - B(k)*temp;
    end
    X(M) = B(M)/D(M);
    for k = M-1:-1:1
        X(k) = (B(k) - X(k+1)*a)/D(k);
    end
    return;

      目前的日期/时间是周一 五月 21, 2018 4:56 pm