MATLAB实例讲解欧拉法求解微分方程

摘要：讲解欧拉法求解微分方程原理，通过MATLAB程序求解实例。

求解微分方程的时候，如果不能将求出结果的表达式，则可以对利用数值积分对微分方程求解，获取数值解。欧拉方法是最简单的一种数值解法。本文理论部分来自知乎作者云端之下的文章“常微分方程——数值解——欧拉方法”，文章链接为：

https://zhuanlan.zhihu.com/p/464118275

实例

求解微分方程dy/dt=-y+t+1,y(0)=1,t的取值为0到2，步长h=0.1,用欧拉法求解微分方程并将结果与y(t)=exp(-t)+t比较。

主程序

clc;
clear all;
close all;
h = 0.1;%步长
y0 = 1;%初值
t = 0:h:2;%x范围
y = exp(-t)+t;%真解
n = length(t);
numy = zeros(1,n);
f2 = -y0;
numy(1) = y0;
%欧拉法计算
for i=2:n
    numy(i) = euler1(y0,h,f2);
    y0 = numy(i);
    f2 = f1(t(i),y0);
end
%绘图
figure;
plot(t,y,'r-','linewidth',1);
hold on;
plot(t,numy,'b-','linewidth',1);
xlabel('t');
grid on;
title('Euler法求系统的输出响应');
ylabel('输出响应y(t)');
legend('真解','Euler法','location','northwest');
wucha_euler = (numy-y).^2;
disp('Euler法误差平方:');
wucha_euler

自定义函数euler1.m

function y = euler1(y0,h,f2)
%%输入参数 y0表示 t=0时 y的取值  即初值
%h表示步长 
%f 表示函数值
%输出y表示方程的响应y
y=y0+h*f2;
end

自定义函数f1.m

function f= f1(x,y)
f = -y+x+1;%微分方程 dy/dt=-y+t+1  初值y(0)=1  微分方程右边的剩余部分构成的函数
end

运行结果

Euler法误差平方:

wucha_euler =

  列 1 至 12

         0    0.0110    0.0097    0.0086    0.0076    0.0067    0.0058    0.0051    0.0044    0.0039    0.0034    0.0029

  列 13 至 21

    0.0025    0.0022    0.0019    0.0016    0.0014    0.0012    0.0010    0.0009    0.0007

改进程序：修改步长h，h分别取值0.1 0.05 0.01 0.001，取值t为0到1，对比结果。

主程序

clc;
clear all;
close all;
h = [0.1 0.05 0.01 0.001];%步长
for j = 1:length(h)
    y0 = 1;%初值
    t = 0:h(j):1;%x范围
    y = exp(-t)+t;%真解
    n = length(t);
    numy = zeros(1,n);
    f2 = -y0;
    numy(1) = y0;
    %欧拉法计算
    for i=2:n
        numy(i) = euler1(y0,h(j),f2);
        y0 = numy(i);
        f2 = f1(t(i),y0);
    end
    %每次因为步长不一样 所以不能用矩阵存结果 
    s{j,:} = numy;%引入元胞类型 Cell 能包含任何类型的数据，比如数值、字符串、逻辑值甚至是Cell自身。
    T{j,:}=t;
end

%绘图
figure;
plot(t,y,'r-','linewidth',1);
hold on;
plot(T{1,:},s{1,:},'b-','linewidth',1);
plot(T{2,:},s{2,:},'g-','linewidth',1);
plot(T{3,:},s{3,:},'k-','linewidth',1);
plot(T{4,:},s{4,:},'m-','linewidth',1);
xlabel('t');
grid on;
title('Euler法求系统的输出响应');
ylabel('输出响应y(t)');
legend('真解','Euler法 h=0.1','Euler法 h=0.05','Euler法 h=0.0.01','Euler法 h=0.001','location','northwest');

自定义函数f1.m

function f= f1(x,y)
f = -y+x+1;%微分方程 dy/dt=-y+t+1  初值y(0)=1  微分方程右边的剩余部分构成的函数
end

自定义函数euler1.m

function y = euler1(y0,h,f2)
%%输入参数 y0表示 t=0时 y的取值  即初值
%h表示步长 
%f 表示函数值
%输出y表示方程的响应y
y=y0+h*f2;
end

运行结果

作 者 | 郭志龙
编 辑 | 郭志龙
校 对 | 郭志龙

本文内容来源于网络，仅供参考学习，如内容、图片有任何版权问题，请联系处理，24小时内删除。