如何用龙格库塔 方法求解 范德波振子的运动
如何用龙格库塔 方法求解 范德波振子的运动
坑，没有想好。
这是一个非常好的线索：
function [t,y] = VanderPol(idis, ivel, mu, A, omega)
[t,y] = ode23(@vdp1,[0 20],[2; 0]);
plot(t,y(:,1),-o,t,y(:,2),-o)
title(Solution of van der Pol Equation (mu = 1) with ODE23);
xlabel(Time t);
ylabel(Solution y);
legend(y_1,y_2)
figure
plot(y(:,1),y(:,2))
title(Phase plane plot)
function dydt = vdp1(t,y)
dydt = [y(2); mu*(1-y(1)^2)*y(2)-y(1)+A*sin(omega*t)];
end
end
在这个代码中，施加了一个外力；
A*sin(omega*t)]
代码这样写，为什们呢？具体有啥规定啊。
---
idis = 2;
ivel = 0;
A = 8;
omega = 4;
mu = 3;
[t,y] = VanderPol(idis, ivel, mu, A, omega);
figure
plot(t,y)
grid
但是这个力不是 岁时间变化的sin函数呢？该怎么办
有想法了：
所有的微分方程组的基本思路：
1、
2、
3、
4、拿 范德波振子举例：
变成了一阶微分方程组之后，
能够非常简单的求得delt x1 和delt x2 如下图
代码如下。
clear all;
x10 = 0.1; %x1代表第一个变量，x2代表第二个变量，后边的0代表的时刻； x20 = 0.1;
n = 1e4; x = zeros(n, 2);
x(1, :) = [x10, x20]; %初始条件；
for k=2:n [dx1, dx2] = dxdt_Vonderpol(x(k-1,1), x(k-1,2)); x(k,1) = x(k-1,1) + dx1; x(k,2) = x(k-1,2) + dx2; end
plot(x(:,1)); grid on;
function g=f1(t, x1, x2) g = x2; end
function g=f2(t, x1, x2) k = 1; %范德波振子的衰减系数 e = 0.086; %其他神经元对当前神经元的耦合项； g = k*(1-x1*x1)*x2 - x1; end
% ============================================================ % --范德波震荡方程的四阶龙格库塔函数 % ============================================================
function [dx1, dx2] = dxdt_Vonderpol(x1, x2) h = 1e-2; %步长
K1 = f1(0, x1, x2); K2 = f1(0, x1 + h*K1/2, x2 + h*K1/2); K3 = f1(0, x1 + h*K2/2, x2 + h*K2/2); K4 = f1(0, x1 + h*K3, x2 + h*K3);
L1 = f2(0, x1, x2); L2 = f2(0, x1 + h*L1/2, x2 + h*L1/2); L3 = f2(0, x1 + h*L2/2, x2 + h*L2/2); L4 = f2(0, x1 + h*L3, x2 + h*L3);
dx1 = (K1 + 2*K2 + 2*K3 + K4)*h/6; dx2 = (L1 + 2*L2 + 2*L3 + L4)*h/6; end
效果如下：
和 微分方程的结果进行对比。
code如下
syms y(t);
[V] = odeToVectorField(diff(y, 2) == (1 - y^2)*diff(y) - y);%把高阶变成一阶；
M = matlabFunction(V,vars, {t,Y});
% sol = ode45(M,[0 30],[2 0]); sol.x 等价于t， sol.y 等价于y； [t, y] = ode45(M,[0 30],[2 0]);
% Plot the Solution y1=y(:,1); %Y y2=y(:,2); %Y
figure(1); plot(t,y1,:b,t,y2,-r) %画微分方程解
figure(2); plot(y1,y2); %画相平面图
结果：
结论，
（1）两种方法的波形是一致的；
（2）初始条件不会影响收敛的相图；
如果加上其他振子的耦合，是类似的，无非 x.. 的微分方程这一块有一点点不一样。
yes
