dy=[y(2);y(3);3*y(3)+y(2)*y(1)]; 注意:尽管不一定用到参数t和y,M一文件必须接受此两参数。这里向量dy必须是列 向量。 (ii)用 Matlab解决此问题的函数形式为 [T, Y]=solver(F', tspan, yO 这里 solver为ode45、ode23、ode113,输入参数F是用M文件定义的常微分方程组, tspan=[ to tfinal]是求解区间,y0是初值列向量。在 Matlab命令窗口输入 [T,Y]=ode45(F',[01],[0;1;-1]) 就得到上述常微分方程的数值解。这里Y和时刻T是一一对应的,Y(:,1)是初值问题的 解,Y(∷,2)是解的导数,Y(∷,3)是解的二阶导数。 例4求 van der pol方程 (1-y2)y+y=0 的数值解,这里>0是一参数。 解(i)化成常微分方程组。设y1=y,y2=y,则有 yI=y y2=(1-y)y2-y1 (i)书写M文件(对于=1) vdl.m function dy=vdp(t, y) dy=[y(2);(1-y(1)^2)*y(2)-y(1)1; (ii)调用Mat1ab函数。对于初值y(0)=2,y(0)=0,解为 [T,Y]=ode45("vdp1',[0201,[2;01) (iv)观察结果。利用图形输出解的结果: title(' Solution of van der Pol Equation, mu=l' xlabel(time t') ylabel('solution y') legend('yl,y2) olution of van der Pol Equation, mu1 例5 van der pol方程,4=1000(刚性) 解(i)书写M文件vdp1000-188- dy=[y(2);y(3);3*y(3)+y(2)*y(1)]; 注意：尽管不一定用到参数t 和 y ，M—文件必须接受此两参数。这里向量 dy 必须是列 向量。 （iii）用 Matlab 解决此问题的函数形式为 [T,Y]=solver('F',tspan,y0) 这里 solver 为 ode45、ode23、ode113，输入参数 F 是用 M 文件定义的常微分方程组， tspan=[t0 tfinal]是求解区间，y0 是初值列向量。在 Matlab 命令窗口输入 [T,Y]=ode45('F',[0 1],[0;1;-1]) 就得到上述常微分方程的数值解。这里 Y 和时刻 T 是一一对应的，Y(:,1)是初值问题的 解，Y(:,2)是解的导数，Y(:,3)是解的二阶导数。 例 4 求 van der Pol 方程 '' (1 ) ' 0 2 y −μ − y y + y = 的数值解，这里 μ > 0是一参数。 解 （i）化成常微分方程组。设 , ' 1 2 y = y y = y ，则有 ⎩ ⎨ ⎧ = − − = 2 1 2 2 1 1 2 ' (1 ) ' y y y y y y μ （ii）书写 M 文件（对于 μ =1）vdp1.m: function dy=vdp1(t,y); dy=[y(2);(1-y(1)^2)*y(2)-y(1)]; （iii）调用 Matlab 函数。对于初值 y(0) = 2, y'(0) = 0 ，解为 [T,Y]=ode45('vdp1',[0 20],[2;0]); （iv）观察结果。利用图形输出解的结果： plot(T,Y(:,1),'-',T,Y(:,2),'--') title('Solution of van der Pol Equation,mu=1'); xlabel('time t'); ylabel('solution y'); legend('y1','y2'); 0 2 4 6 8 10 12 14 16 18 20 -3 -2 -1 0 1 2 3 Solution of van der Pol Equation,mu=1 time t solution y y1 y2 例 5 van der Pol 方程， μ =1000 （刚性） 解 (i)书写 M 文件 vdp1000.m: