考虑如下的一阶微分方程:
$\frac{dx}{dt}=f(x,t)$
选取小的h,有如下结果:
$\frac{dx}{dt}=f(x,t)=\frac{x(t+h)-x(t)}{h}=\frac{x_k-x_{k-1}}{h}$
重写为如下形式:
$x_k=x_{k-1}+hf(x,t)$
使用下面的二阶线性微分方程作为例子:
$\frac{d^2x}{dt^2}=-\omega ^2x$
初始条件:
$x(0)=0$
$\frac{dx}{dt}=\omega$
其解析解为$x=\sin{\omega t}$
matlab程序如下(请忽视作为一个初学者的糟糕语法)
% Euler numerical integration for d(dx/dt)/dt = -omega*omega*x
w = 2; % w is omega
x = zeros(1, 1000); % x init condition:=0
xd = w; % xd = dx/dt, init condition:=w
h=0.01; % step size
loop_times = 1;
for t = 0:h:9.99
xdd = -w * w * x(loop_times);
xd = xd + h*xdd;
x(loop_times + 1) = x(loop_times) + h*xd;
loop_times = loop_times + 1;
end
plot(x);
hold on; % put the following curve on the same plot
t = 0:h:9.99;
plot(sin(w*t));