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

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

实例

求解微分方程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表示方程的响应yy=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表示方程的响应yy=y0+h*f2;end

运行结果

作 者 | 郭志龙

编 辑 | 郭志龙

校 对 | 郭志龙

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