求解微分方程的时候，如果不能将求出结果的表达式，则可以对利用数值积分对微分方程求解，获取数值解。欧拉方法是最简单的一种数值解法。前面介绍过MATLAB实例讲解欧拉法求解微分方程，今天实例讲解欧拉法求解一阶微分方程组。

本文理论部分来自知乎作者云端之下的文章“常微分方程——数值解——欧拉方法”，文章链接为：

实例

假设一个一阶的动态系统为：

dx/dt = y

dy/dt = 0.04y-1.2x

dz/dt = 0.25

动态系统的初始条件t = 0时刻为：x(0) = 2.0，y(0) = 0.0，z(0) = 0.0。使用欧拉法，得到t从0到15时刻所有x、y和z的值，步长是h =0.02s。

主程序

clc;clear all;close all;h = 0.02;%步长y0 = [2 0 0];%初值t = 0:h:15;%x范围n = length(t);numy = zeros(n,3);f2 = y0;numy(1,:) = y0;%欧拉法计算for i=2:n    numy(i,:) = euler11(y0,h,f2);    y0 = numy(i,:);    f2 = f11(t(i),y0);endx = numy(:,1);y = numy(:,2);z = numy(:,3);figure;plot3(x,y,z,'k--','linewidth',1);xlabel('x');ylabel('y');zlabel('z');title('solution');grid on;

euler11.m自定义函数

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

f11.m自定义函数

function f= f11(t,Y)% y(1):x y(2):y  y(3):zx = Y(1);y = Y(2);z = Y(3);f(1) = y;f(2) = 0.04*y-1.2*x;f(3) = 0.25;end

运行结果

方式二程序

主程序

clc;clear all;close all;h = 0.02;%步长y0 = [2 0 0];%初值t = 0:h:15;%x范围n = length(t);numy = zeros(n,3);f2 = y0;numy(1,:) = y0;%欧拉法计算for i=2:n    numy(i,:) = euler11(y0,h,f2);    y0 = numy(i,:);%     f2 = f11(t(i),y0);    f2 = f12(t(i),y0);endx = numy(:,1);y = numy(:,2);z = numy(:,3);figure;plot3(x,y,z,'k--','linewidth',1);xlabel('x');ylabel('y');zlabel('z');title('solution');grid on;

euler11.m自定义函数

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

f12.m程序

function f = f12(t,y)% y(1):x y(2):y  y(3):zf(1) = y(2);%微分方程右边的剩余部分构成的函数f(2) = 0.04*y(2)-1.2*y(1);f(3) = 0.25;end

运行结果

作 者 | 郭志龙

编 辑 | 郭志龙

校 对 | 郭志龙

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