MATLAB在求解微分方程的时候怎么同时输出所求函数的导数？（使用ode45的同时输出微分方程结果）

以Van Del Pol这样一个微分方程为例：

求解这样一个微分方程，MATLAB代码如下：

tspan=0:0.01:15;
q0=[2;0];%初值 [q;q']
[t,q]=ode45(@vdp1,tspan,q0);

plot(t,q(:,1),t,q(:,2))
title('Solution of van der Pol Equation (\mu = 1) with ODE45');
xlabel('Time t');
ylabel('Solution q dq');
legend('q_1','q_2')


%将以下代码保存为vdp1.m
function qdot=vdp1(t,q)
%求Van Der Pol 方程
% 令q1=q' ,q2=q'' 把q''+(q^2-1)q'+q=0写成状态方程
qdot=zeros(2,1);%使xdot成为二元列向量
qdot(1)=q(2);
qdot(2)=(1-q(1)^2)*q(2)-q(1);

end

结果如图：

所获得的结果是q 和 dq

这时候我们想获得ddq的结果，该怎么办呢？

有两个办法：
1.使用deval函数
2.使用gradient函数

下面分别叙述

1.使用deval函数

tspan=0:0.01:15;
q0=[2;0];%初值 [q;q']
sol=ode45(@vdp1,tspan,q0);
[q,dq]=deval(sol,tspan)
plot(t,q(1,:),t,q(2,:),'ro')
title('Solution of van der Pol Equation (\mu = 1) with ODE45');
xlabel('Time t');
ylabel('Solution q');
legend('q','q''')

figure(2)
plot(t,dq(1,:),'ro',t,dq(2,:))
xlabel('Time t');
ylabel('Solution dq');
legend('q''','q''''')


%以下代码保存为vdp1.m
function qdot=vdp1(t,q)
%求Van Der Pol 方程
% 令q1=q' ,q2=q'' 把q''+(q^2-1)q'+q=0写成状态方程
qdot=zeros(2,1);%使xdot成为二元列向量
qdot(1)=q(2);
qdot(2)=(1-q(1)^2)*q(2)-q(1);

end

结果如图

输出的结果是 dq和ddq
从dq结果的对比可以看到，此方法有效。

2.使用gradient函数

tspan=0:0.01:15;
q0=[2;0];%初值 [q;q']
[t,q]=ode45(@vdp1,tspan,q0);
[~,dq] = gradient(q, mean(diff(t)));

plot(t,q(:,1),t,q(:,2),'ro')
title('Solution of van der Pol Equation (\mu = 1) with ODE45');
xlabel('Time t');
ylabel('Solution q');
legend('q','q''')

figure(2)
plot(t,dq(:,1),'ro',t,dq(:,2))
xlabel('Time t');
ylabel('Solution dq');
legend('q''','q''''')



%以下代码保存为vdp1.m
function qdot=vdp1(t,q)
%求Van Der Pol 方程
% 令q1=q' ,q2=q'' 把q''+(q^2-1)q'+q=0写成状态方程
qdot=zeros(2,1);%使xdot成为二元列向量
qdot(1)=q(2);
qdot(2)=(1-q(1)^2)*q(2)-q(1);

end

输出结果如图所示：

同样，gradient也是实现输出所要求解微分方程的导数。

从下图可以看到:两个函数实现的效果是一样的.