摘要：讲解欧拉法求解微分方程原理，通过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,&39;r-&39;,&39;linewidth&39;,1);hold on;plot(t,numy,&39;b-&39;,&39;linewidth&39;,1);xlabel(&39;t&39;);grid on;title(&39;Euler法求系统的输出响应&39;);ylabel(&39;输出响应y(t)&39;);legend(&39;真解&39;,&39;Euler法&39;,&39;location&39;,&39;northwest&39;);wucha_euler = (numy-y).^2;disp(&39;Euler法误差平方:&39;);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,&39;r-&39;,&39;linewidth&39;,1);hold on;plot(T{1,:},s{1,:},&39;b-&39;,&39;linewidth&39;,1);plot(T{2,:},s{2,:},&39;g-&39;,&39;linewidth&39;,1);plot(T{3,:},s{3,:},&39;k-&39;,&39;linewidth&39;,1);plot(T{4,:},s{4,:},&39;m-&39;,&39;linewidth&39;,1);xlabel(&39;t&39;);grid on;title(&39;Euler法求系统的输出响应&39;);ylabel(&39;输出响应y(t)&39;);legend(&39;真解&39;,&39;Euler法 h=0.1&39;,&39;Euler法 h=0.05&39;,&39;Euler法 h=0.0.01&39;,&39;Euler法 h=0.001&39;,&39;location&39;,&39;northwest&39;);

自定义函数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小时内删除。