[数学]10实验十 常微分方程.ppt

【例6】求初值问题 在区间[1.2,4]上的近似解, 并作图。 输入： fun=inline((1+x*y)*y/(x*y-1),x,y); [x,y]=ode45(fun,[1.2,4],1); [x,y] ode45(fun,[1.2,4],1) %或plot(x,y) 其输出为数值近似解(插值函数)的形式： ans = 1.2000 1.0000 1.2046 1.0456 1.2091 1.0845 1.2137 1.1194 1.2183 1.1518 1.2274 1.2114 1.2366 1.2658 1.2457 1.3168 1.2549 1.3654 …… 3.8290 35.1151 3.8860 37.2059 3.9430 39.4189 4.0000 41.7613 并得到近似解的图形(见图10.5)。 图10.5 【例7】求VanderPol方程 , , 在区间[0,20]上的近似解。 令 , , 原方程化为方程组。 单击File：New：M-file打开命令窗口, 输入以下程序, 建立并保存函数文件vdp1.m： function dydt=vdp1(t,y) dydt=[y(2);(1-y(1)^2)*y(2)-y(1)]; 回到MTLAB命令窗口, 输入求VanderPol方程近似解命令： [t,y]=ode23(@vdp1,[0 20],[0 -0.5]); plot(t,y(:,1)); 可以观察到近似解的图形(见图10.6)。 图10.6 10.2.4 一阶微分方程的斜率场 设一阶微分方程的形式为 , 即研究导数已经解出的一阶方程。其中 是已知函数。由导数的几何意义, 未知函数y的斜率是平面上点(x,y)的函数。因此可以在xoy平面上的每一点, 作出过该点的以 为斜率的一条很短的直线(即是未知函数y的切线)。这样得到的一个图形就是微分方程 的斜率场。为了便于观察, 实际上只要在xoy平面上取适当多的点, 作出在这些点的函数的切线。观察这样的斜率场, 就可以得到方程 的近似积分曲线。 【例8】绘制微分方程 的斜率场和解曲线。 因为没有给出初始条件, 所以先建立并保存成函数文件zxyf.m。 function Dy=zxyf(x,y) Dy=x./y; %计算斜率函数 编写命令文件： clf, clear; a=0;b=4;c=0;d=4; n=15; [x,y]=meshgrid(linspace(a,b,n),linspace(c,d,n)); %生成区域中的网格 f=x./y; Fx=cos(atan(x./y)); Fy=sqrt(1-Fx.^2); %计算切线斜率矢量 quiver(x,y,Fx,Fy,0.5) %在每一网格点画出相应的切线斜率矢量,0.5是控制矢量大小的参数 hold on axis([a,b,c,d]) 输出见图10.7。 图10.7 如果增加初始条件y(0)=2, 则得到的是一条确定的积分曲线。输入： [x,y]=ode45(zxyf,[0,4],2); plot(x,y,r.-) 其输出为图10.8中的曲线(屏幕上显示为红色)。 图10.8 【例9】输入以下命令, 观察所得到的斜率场和等高线： clf [x,y]=meshgrid(-2:0.2:0,-1:0.15:1); z=x.^exp(-x.^2+y.^2); [px,py]=gradient(z,0.2,0.15); contour(x,y,z) %显示矩阵z的二维等高线图 hold on quiver(x,y,px,py) hold off axis image 输出见图10.9。 图10.9 MATLAB 高等数学实验 实验十 常微分方程 实验目的 掌握用MATLAB求微分方程及方程组解析解的方法, 学习求微分方程近似解的欧拉折线法和龙格-库塔法。 10.1 学习MATLAB命令 10.1.1 求常微分方程的符号解函数 dsolve(eq1,eq2,…, cond1,cond2,…, v) 输入符号形式的常微分方程(组) eq1