梯形方法求解常微分方程

Mi@MI

已于 2023-11-24 21:31:12 修改

阅读量1k

点赞数 12

分类专栏：计算数学|数值解文章标签：算法线性代数数学建模

于 2023-11-24 21:30:11 首次发布

本文链接：https://blog.csdn.net/m0_67196062/article/details/134605195

版权

计算数学|数值解专栏收录该内容

7 篇文章 1 订阅

订阅专栏

对于一阶常微分方程初值问题：

$\left\{\begin{matrix} \frac{dy}{dx}=f(x,y),a<x<b \\ y(a)=\alpha \end{matrix}\right.$

利用数值积分可以推导出如下的梯形方法：

$y_{i+1}=y_{i}+h/2[f(x_{i},y_{i})+f(x_{i+1},y_{i+1})]$

从公式可以看出，梯形方法是隐式方法，因此需要迭代法进行求解，这里迭代初值采用欧拉方法，当然也可以采用其他方法，只要是显式方法即可，采用欧拉方法提供初值的改进梯形公式如下：

$\left\{\begin{matrix} y_{i+1}^{(0)}=y_{i}+hf(x_{i},y_{i})\\ y_{i+1}^{(k+1)}=y_{i}+h/2[f(x_{i},y_{i})+f(x_{i+1},y_{k}^{(k)})] \end{matrix}\right.$

使用迭代法需要设置误差限以控制迭代步数，数值算例如下：

$\left\{\begin{matrix} \frac{dy}{dx}=y-\frac{2x}{y},0< x\leq 1\\ y(0)=1 \end{matrix}\right.$

步长取0.1，迭代误差限取 $\varepsilon =10^{-4}$ ，MATLAB主函数代码如下：

% 梯形方法，为隐式方法，需要借助欧拉方法进行迭代
% L区间右端点，左端点默认为0，h为步长
function [err]=tixing(L,h)
x=(0:h:L);
% 设置误差限
eps=10^(-4);
%初始化分配内存
ye=zeros(1,length(x));y=ye;
%初值
y(1)=1;
% 精确解
for i=1:length(x)
    ye(i)=sqrt(1+2*x(i));
end
% 数值解
% 先用欧拉方法提供初值,再利用梯形方法
for i=1:length(x)-1
y1=y(i)+h*(y(i)-2*(x(i)/y(i)));
y2=y(i)+h/2*(y(i)-2*(x(i)/y(i))+y1-2*(x(i+1)/y1));
    while 1
        if abs(y1-y2)<eps
            y(i+1)=y2;
            break
        else
            y1=y2;
            y2=y(i)+h/2*(y(i)-2*(x(i)/y(i))+y1-2*(x(i+1)/y1));
            continue
        end

    end
end
% 计算误差
err=abs(ye-y);

计算收敛阶代码如下：

% 梯形方法收敛阶计算
dh=[1/100,1/200,1/400,1/800];
n=length(dh);
for i=1:n
    e(i)=max(tixing(1,dh(i)));
end
for i=1:n-1
    rate_tixing(i)=log2(e(i)/e(i+1));
end
% 由结果可以知道梯形方法是二阶收敛的

并且结合主函数代码可以得出梯形方法确实二阶收敛的，与理论分析是一致的，另外考虑到也有些学者们使用python编程，因此这里也给出python代码，如下：

import numpy as np
def euler_tixing(f,x,y0,args=(),eps=1e-4,max_iter=100):
    """
    :param f: 右端函数
    :param x: 离散的节点
    :param y0: 初值
    :param args:
    :param eps: 精度
    :param max_iter: 最大迭代次数
    :return: 返回数值解
    """
    y=np.zeros_like(x)
    y[0]=y0
    for i in range(len(x)-1):
        h=0.1
        y[i+1]=y[i]+h*f(x[i],y[i],*args)
        for j in range(max_iter):
            last_y=y[i+1]
            y[i+1]=y[i]+0.5*h*(f(x[i],y[i],*args)+f(x[i+1],y[i+1],*args))
            if np.abs(y[i+1]-last_y)<eps:
                break
            else:
                pass
    return y
x=np.linspace(0,1,num=11)
y0=1
y=euler_tixing(lambda x,y:y-2*x/y,x,y0)
print('   步长，   数值解，    精确解，     误差')
for i in range(len(x)):
    ans=np.sqrt(1+2*x[i])
    err=np.abs(y[i]-ans)
    print(f'{x[i]:4.0f},{y[i]:10.6f},{ans:10.6f},{err:10.2e}')

至于python计算收敛阶，这里就留作一个小习题，感兴趣的小伙伴可以尝试尝试

PS：最后，初来乍到，觉得满意的小伙伴们能否点个赞和关注呢，你们的支持将是我创作的最大动力