数字积分

考虑如下的一阶微分方程:

$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));

原文地址:https://www.cnblogs.com/byeyear/p/6714342.html