数值解求解函数有[t,x].pptVIP

ODE的求解 Contents Contents 1.常微分方程概念 常微分方程概念（续） 常微分方程概念（续） 常微分方程概念（续） 2.常微分方程的应用 例.几个常见的微分方程模型 续上例 结果很复杂，对其进行化简： x=simple(x) y=simple(y) 得到结果： x=5/7*exp(-t)-9/48*C1*exp(11/12*t+1/12*t*193^(1/2)) +1/48*C1*exp(11/12*t+1/12*t*193^(1/2))*193^(1/2)-49/48*C2*exp(11/12*t-1/12*t*193^(1/2))-1/48*C2*exp(11/12*t-1/12*t*193^(1/2))*193^(1/2) y=-8/7*exp(-t)+C1*exp(11/12*t+1/12*t*193^(1/2))+ C2*exp(11/12*t-1/12*t*193^(1/2)) 还可以得到近似的结果： x=vpa(x，4) %保留4位有效数字 y=vpa(y，4) x =.7143*exp(-1.*t)-.7314*C1*exp(2.074*t)-1.310*C2*exp(-.2410*t) y =-1.143*exp(-1.*t)+.9998*C1*exp(2.074*t)+.9998*C2*exp(-.2410*t) 其中，C1、C2为Matlab给出的任意常数. 例：高等数学中的例题 例2.P278 一电路（电源串联一电阻R和一电感L），电源电动势为E=Emsin wt，电阻R和电感L为常量，求电流I（t）. 解：显然电流满足如下常微分方程： \omega用w代替 DI+R/L*I=Em/L*sin wt 设电路接通时刻为t0=0，于是有初始条件I(0)=0 在Matlab下求解 syms R L Em w; dsolve(‘DI+R/L*I=Em/L*sin(w*t)’,‘I(0)=0’)； 得到 ans =exp(-R/L*t)/(R^2+w^2*L^2)*Em*w*L-Em*(w*L*cos(w*t)-R*sin(w*t))/(R^2+w^2*L^2) 键入pretty(simple(ans))得到 例：高等数学中的例题 P318 求Euler方程的通解 求解命令 dsolve(x^3*D3y+x^2*D2y-4*x*Dy=3*x^2,x) 得到通解 1/3*C1*x^3-1/2*x^2-C2/x+C3 与用变换和特征方程的方法求解相比…… 例.数学摆的求解 在数学摆模型中考虑m=1kg,l=1m,g=9.8m/s^2,初始角度\theta_0=pi/12，初速度为0，阻尼系数\lambda=0.1 根据前面描述无阻尼数学摆的方程可以化为如下方程组： 可以得到数值解图像 键入：theta0=[pi/12,0];[t,theta]=ode23(fun2,[0,100],theta0); plot(t,theta(:,1)); 例：Voterra系统 %……………………vot.m…………………… function xdot=vot(t,x) a=1;c=-.1;d=-.5;e=.02;f=0;b=0; xdot=[x(1)*(a+b*x(1)+c*x(2));x(2)*(d+e*x(1)+f*x(2))]; %……………………vot.m…………………… [t,x]=ode45(vot,[0,300],x0) plot(x(:,1),x(:,2)) 例：Logistic 模型决定的向量场 （续上例）参考程序 clear c=0.01; x0=0:4:100; y0=0:4:100; [x,y]=meshgrid(x0,y0); %设定绘制向量场的点 d=sqrt(1+(0.1*(1-y/50).*y).^2); u=c./d; %与v一同决定向量场方向 v=c*(0.1*(1-y/50).*y)./d; hold on quiver(x,y,u,v); 例.简单的差分方程 模拟斐波那契数列(第一个月1对小兔子开始，每个月的兔子数量) %rabbit.mdl 例 例 求解Dy=-\alpha*y,y(0)=100 例.求数学摆的数值解（无阻尼） 例.求数学摆的数值解（有阻尼） 求解*：Dx=-0.1*t*x 作业II 1.给出适当参数绘制出SIR模型的解并给出适当的解释; 2.#熟悉ODE的解析解与数值解的Matlab求解方法； 3.了解一维方程向量场绘图；自己给出适当方程，绘图、观察积分曲线的走向； 4.（1）求解Van de pol方程，当$\mu=1$时； （2）*尝试在matlab下输入以下命令： [t,y]=ode45(@vdp1,[0 20],