以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也是实现输出所要求解微分方程的导数。
从下图可以看到:两个函数实现的效果是一样的.