实验二微分方程与差分方程模型matlab求解.docVIP

实验二： 微分方程与差分方程模型Matlab求解 一、实验目的[1] 掌握解析、数值解法，并学会用图形观察解的形态和进行解的定性分析；[2] 熟悉MATLAB软件关于微分方程求解的各种命令； [3] 通过范例学习建立微分方程方面的数学模型以及求解全过程； [4]Logistic模型的求解与混沌的产生过程。 二、实验原理1 求解一阶微分方程 ??? (1) 求通解 输入： dsolve(Dy=1+y^2) 输出： ans = tan(t+C1) （2）求特解 输入： dsolve(Dy=1+y^2,y(0)=1,x) 指定初值为1，自变量为x 输出： ans = tan(x+1/4*pi) 例2 求解二阶微分方程 原方程两边都除以，得 输入: dsolve(D2y+(1/x)*Dy+(1-1/4/x^2)*y=0,y(pi/2)=2,Dy(pi/2)=-2/pi,x) ans = - (exp(x*i)*(pi/2)^(1/2)*i)/x^(1/2) + (exp(x*i)*exp(-x*2*i)*(pi/2)^(3/2)*2*i)/(pi*x^(1/2)) 试试能不用用simplify函数化简 输入: simplify(ans) ans = 2^(1/2)*pi^(1/2)/x^(1/2)*sin(x) ? （2）微分方程组 例3 求解???? df/dx=3f+4g;? dg/dx=-4f+3g。 f = exp(3*t)*(C1*sin(4*t)+C2*cos(4*t)) g = exp(3*t)*(C1*cos(4*t)-C2*sin(4*t)) 特解： [f,g]=dsolve(Df=3*f+4*g,Dg=-4*f+3*g,f(0)=0,g(0)=1) f = exp(3*t)*sin(4*t) g = exp(3*t)*cos(4*t) 数值解 在微分方程(组)难以获得解析解的情况下，可以用Matlab方便地求出数值解。格式为: [t,y] = ode23(F,ts,y0,options) 注意： 微分方程的形式：y = F(t, y)，tts=[t0, tf] 表示自变量的取值范围，（2）ts=[t0,t1,t2,…,tf],则输出在指定时刻t0,t1,t2,…,tf处给出，（3）ts=t0:k:tf,则输出在区间[t0,tf]的等分点给出； y0为初值条件； options用于设定误差限（缺省是设定相对误差是10^(-3)，绝对误差是10^(-6)） 解 形式转化：令。则以上方程转化一阶微分方程组： 。 编写M文件如下，必须是M文件表示微分方程组，并保存，一般地，M文件的名字与函数名相同，保存位置可以为默认的work子目录，也可以保存在自定义文件夹，这时注意要增加搜索路径（File\Set Path\Add Folder） ??? dot1=[y(2); (1-y(1)^2)*y(2)-y(1)]; 在命令窗口写如下命令: [t,y]=ode23(vdpol,[0,20],[1,0]); y1=y(:,1);y2=y(:,2); plot(t,y1,t,y2,--);title(Van Der Pol Solution ); xlabel(Time,Second);ylabel(y(1)andy(2)) 执行： 注：Van der Pol方程描述具有一个非线性振动项的振动子的运动过程。最初,由于它在非线性电路上的应用而引起广泛兴趣。一般形式为。 图形解 无论是解析解还是数值解，都不如图形解直观明了。即使是在得到了解析解或数值解的情况下，作出解的图形，仍然是一件深受欢迎的事。这些都可以用Matlab方便地进行。 (1)图示解析解 如果微分方程（组）的解析解为：y=f (x)，则可以用Matlab函数fplot作出其图形： fplot(fun,lims) 其中：fun给出函数表达式；lims=[xmin xmax ymin ymax]限定坐标轴的大小。例如 fplot(sin(1/x), [0.01 0.1 -1 1]) (2)图示数值解 设想已经得到微分方程（组）的数值解(x,y)。可以用Matlab函数plot(x,y)直接作出图形。其中x和y为向量（或矩阵）。 2、Volterra模型（食饵捕食者模型） 意大利生物学家Ancona曾致力于鱼类种群相互制约关系的研究，他从第一次世界大战期间,地中海各港口捕获的几种鱼类捕获量百分比的资料中，发现鲨鱼的比例有明显增加（见下表）。 年代 1914 1915 1916 1917 1918 百分比 11.9 21.4 22.1 21.2 36.4 年代 1919 1920 1921 1922 1923