计算物理学（数值分析）上机实验答案5、常微分方程初值问题的数值解法

实验五、常微分方程初值问题的数值解法

​ 常微分方程的求解问题在实践中经常遇到，因此研究常微分方程的数值解法就很有必要。欧拉方法是最简单、最基本的方法，利用差商代替微商，就可得到一系列欧拉公式。这些公式形式简洁，易于编程计算，但精度较低，可方便用于精度要求不高的近似计算。龙格-库塔方法是利用区间上多个点的斜率值的加权平均的思想，得出的高精度的计算公式。特别是四阶龙格-库塔公式，易于编程计算，精度较高，是常用的工程计算方法。线性多步方法是在用插值多项式代替被积函数的基础上构造的，它可利用前面若干步计算结果的信息，使计算结果提高精度，其中较常用的是四阶阿当姆斯方法。

一、实验目的

1、 熟悉求解常微分方程初值问题的有关方法和理论，主要是改进欧拉法、四阶龙格-库塔法；

2、会编制上述方法的计算程序，包括求解微分方程组的计算程序；

3、针对实习题编制程序，并上机计算其所需要的结果；

4、 通过对各种求解方法的计算实习，体会各种解法的功能、优缺点及适用场合， 会选取适当的求解方法。

二、算法实例

1.向前欧拉公式的 MATLAB 主程序

例 5.1

​ 用欧拉方法求初值问题的数值解，分别取h = 0.0750，0.0075，并计算误差，画出精确解和数值解的图形。
$$\{\begin{array}{c}\frac{dy}{dx}=x-y,0\le x\le 1\\ y{|}_{x=0}=1\end{array}$$

解：编写并保存名为 Eulerli1.m 的 MATLAB 计算和画图的主程序如下

function P=Eulerli1(x0,y0,b,h)
n=round(fix((b-x0)/h));%取整数
X=zeros(n,1);Y=zeros(n,1); k=1;
X(k)=x0; Y(k)=y0;
for k=1:n
    X(k+1)=X(k)+h;
    Y(k+1)=Y(k)+h*(X(k)-Y(k)); k=k+1;
end
y=X-1+2*exp(-X); plot(X,Y,'mp',X,y,'b-')
grid
xlabel('自变量 X'), ylabel('因变量 Y')
title('用向前欧拉公式求dy/dx=x-y，y(0)=1在[0,1]上的数值解和精确解y=x-1+2 exp(-x)')
legend('h=0.075时，dy/dx=x-y，y(0)=1在[0,1]上的数值解','精确解y=x-1+2 exp(-x)')
jwY=y-Y;xwY=jwY./y;
k1=1:n;k=[0,k1];
P=[k',X,Y,y,jwY,xwY];

在 MATLAB 工作窗口输入下面的程序

x0=0;y0=1;b=1;h=0.0750;
P=Eulerli1(x0,y0,b,h)

在 MATLAB 工作窗口输入下面的程序

h1=0.0075; P1=Eulerli1(x0,y0,b,h1)
legend('h1=0.0075 时，dy/dx=x-y，y(0)=1 在[0,1]上的数值解','精确解 y=x-1+2 exp(-x)')

2.自适应向前欧拉公式的 MATLAB 主程序

function [H,X,Y,k,h,P]=QEuler(funfcn,x0,b,y0,tol)
%初始化.
pow=1/3;
if nargin < 5 | isempty(tol)
    tol = 1e-6;
end;
if nargin < 6 | isempty(trace)
    trace = 0;
end;
x=x0; h=0.0078125*(b-x); y=y0(:);p=128;
H=zeros(p,1);
X=zeros(p,1); Y=zeros(p,length(y)); k=1; X(k)=x;
Y(k,:)=y';
% 绘图.
clc,x,h,y
end
% 主循环.
while (x<b)&(x+h>x)
    if x+h>b
        h=b-x;
    end
    %计算斜率.
    fxy=feval(funfcn,x,y); fxy=fxy(:);
    %计算误差，设定可接受误差.
    delta=norm(h*fxy,'inf');
    wucha=tol*max(norm(y,'inf'),1.0);
    % 当误差可接受时重写解.
    if delta<=wucha
        x=x+h; y=y+h*fxy; k=k+1;
        if k>length(X)
            X=[X;zeros(p,1)]; Y=[Y;zeros(p,length(y))];
            H=[H;zeros(p,1)];
        end
        H(k)=h;X(k)=x;Y(k,:)=y'; plot(X,Y,'rp'), grid
        xlabel('自变量 X'), ylabel('因变量 Y')
        title('用向前欧拉（Euler）公式计算dy/dx=f(x,y)，y(x0)=y0在[x0,b]上的数值解')
    end
    %更新步长.
    if delta~=0.0
        h=min(h*8,0.9*h*(wucha/delta)^pow);
    end
end
if (x<b)
    disp('Singularity likely.'), x
end
H=H(1:k);
X=X(1:k);
Y=Y(1:k,:);
n=1:k; P=[n',H,X,Y]

3.向后欧拉方法的 MATLAB 程序

function [X,Y,n,P]=Heuler1(funfcn,x0,b,y0,h,tol)
n=fix((b-x0)/h); X=zeros(n+1,1); Y=zeros(n+1,1);
k=1; X(k)=x0; Y(k,:)=y0; Y1(k,:)=y0;
% 绘图.
clc,x0,h,y0
% 产生初值.
for i=2:n+1
    X(i)=x0+h; Y(i,:)=y0+h*feval(funfcn,x0,y0);
    Y1(i,:)=y0+h*feval(funfcn,X(i),Y(i,:));
    % 主循环.
    Wu=abs(Y1(i,:)-Y(i,:));
    while Wu>tol
        p=Y1(i,:);
        Y1(i,:)=y0+h*feval(funfcn,X(i),p);
        Y(i,:)=p;
    end
    x0=x0+h; y0=Y1(i,:);
    Y(i,:)=y0; plot(X,Y,'ro')
    grid on
    xlabel('自变量 X'), ylabel('因变量 Y')
    title('用向后欧拉公式计算dy/dx=f(x,y)，y(x0)=y0在[x0,b]上的数值解')
end
X=X(1:n+1); Y=Y(1:n+1,:); n=1:n+1; P=[n',X,Y]

例 5.2

​ 用向后欧拉公式求解区间[0,2] 上的初值问题$\frac{dy}{dx}=-3y+8x-7$