数值计算讨论组

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


    [微]显式Euler方法,梯形方法,预估-校正Euler方法参考程序及解释

    分享

    Admin
    Admin

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

    [微]显式Euler方法,梯形方法,预估-校正Euler方法参考程序及解释

    帖子 由 Admin 于 周二 十月 26, 2010 10:40 pm

    程序中的注释充分详细。大部分同学没有解决好 梯形方法中求隐式方程的部分,请参考程序中的while循环。改程序给出了5种不同步长下的求解区间右端点处的误差值,分别存在了errF,errT,errPC中,程序运行之后会给出errF,errT,errPC及5幅对比图像。通过errF中误差随步长减半而减半可知显式Euler方法整体截断误差阶为1;errT和errPC中的数值表明,步长减半时误差逐减为原1/4,可知梯形方法和预估校正Euler方法的整体截断误差阶为2。


    代码:
    % 本程序实现求解下述问题的 显式Euler方法,梯形Euler方法,预估校正Euler方法
    % 一阶常微分方程初值问题: y' = - y^2,  y(0)=1, [0,1], 精确解y(x)=1/(1+x)

    clear;clc;close all;
    format long
    % errF,errT,errPC分别记录三种方法在5种步长下 求解区间右端点处的误差
    errF = zeros(5,1);
    errT = zeros(5,1);
    errPC = zeros(5,1);
    n = 8;
    for M = 1:5
        % 用下式来给出不同的步长h=1/6,1/32,1/64,1/128,1/256
        n = n*2;
        h = 1/n;
        x = 0:h:1;
           
        % 精确解
        y_exact = zeros(n+1,1);
        for i = 1:n+1
            y_exact(i) = 1/(1+x(i));
        end

        % 显式Euler方法
        yF = zeros(n+1,1);
        yF(1) = 1;
        for i = 1:n
            yF(i+1) = yF(i) - h*yF(i)^2;
        end

        % 梯形方法
        yT = zeros(n+1,1);
        yT(1) = 1;
        for i = 1:n
            % 以这种方式给初值
            yT(i+1) = yT(i) - h*yT(i)^2;
            % 简单迭代求解隐式方程
            err = 1;
            k = 0;
            while err> 1.0e-9
                temp = yT(i) - h*(yT(i+1)^2 + yT(i)^2)/2;
                err = abs(temp - yT(i+1));
                yT(i+1) = temp;
                % 用k记录迭代次数
                k = k + 1;
            end
        end

        % 预估-校正Euler方法
        yPC = zeros(n+1,1);
        yPC(1) = 1;
        for i = 1:n
            yPC(i+1) = yPC(i) - h*yPC(i)^2;
            yPC(i+1) = yPC(i) - h*(yPC(i+1)^2 + yPC(i)^2)/2;
        end

        figure
        plot(x,y_exact,'r',x,yF,'-.',x,yT,'+',x,yPC,'*')
        legend('精确','显式Euler','梯形','预校Euler')
       
       
        errF(M) = abs(yF(n+1) - y_exact(n+1));
        errT(M) = abs(yT(n+1) - y_exact(n+1));
        errPC(M) = abs(yPC(n+1) - y_exact(n+1));
    end
    errF
    errT
    errPC


      目前的日期/时间是周二 七月 17, 2018 8:01 am