2016春 作业 实验1常微分方程
常微分方程数值解实验
有大量的常微分方程虽然从理论上讲，其解是存在的，但我们却无 法求出其解析解，此时，我们需要寻求方程的数值解，在求常微分方程 数值解方面，MATLAB具有丰富的函数，我们将其统称为solver，其一般 格式为：
Image
Image
如果微 分方程 由一个 或多个 高阶常微分方程给出，要得到该方程的数值解,可以将方程转换成一阶 常微分方程组。假设高阶常微分方程的一般形式为y( n) = f ( t, y, yʹ, ⋯,y( n - 1) )，而且函数y(t)的各阶导数初值为y(0)，yʹ(0) ,…, y( n - 1) (0)可以选 择一组变量令: x1= y, x2 = yʹ,…, xn = y( n - 1) ,我们就可以把原高阶常微 分方程转换成下面的一阶常微分方程组形式： 而且初值x1(0)=y(0),x2(0)=yʹ(0),…，xn(0)=(0)。 转换以后就可以求原 高阶常微分方程的数值解了。 例2 求微分方程，,的数值解。 对方程进行变换,选择变量 (1) 建立自定义函数 用edit命令建立自定义函数名为f.m,内容为: function y =f(t,x) y=[x(2);x(3);-t^2*x(2)*x(1)^2-t*x(1)*x(3)+exp(t*x(1))]; (2)调用对微分方程数值解ode45函数求解 用edit命令建立一个命令文件f2. m,内容为: >>x0=[2;0;0]; >>[t,y] =ode45(’f’,[0,10],x0);plot(t,y); >>figure; >>plot3(y(:,1),y(:,2), y(:,3))得
常微分方程作业
常微分⽅程作业安顺市镇宁县六马中学教师：韦应俭第⼀部分⼀、常微分⽅程的概念含有⾃变量、函数及其导数的关系式. ⼆、⼀阶微分⽅程的初等解法（1）变量分离⽅程形如：)()(y x f dydxρ=的⽅程，称为变量分离⽅程，这⾥)(),(y x f ρ分别是y x ,的连续函数.（2）可化为分离变量⽅程的⽅程的三种形式①)(xy f dy dx yx =?；②)(x y g dy dx =；③)(222111xc x b x a x c x b x a f dy dx++++= （3）贝努⼒⽅程n y x g y x dydx)()(+=ρ（4）⼀阶线性⽅程)()(x g y x dxdy+=ρ（5）Riccaiti ⽅程)()()(2x r y x g y x dxdy++=ρ（6）形如0),(),(=+dy y x N dx y x M 的⽅程①若0=??-??xNy M ，则⽅程式恰当的通解是0)(.0)1(12=-+==+-+dy x y ydx dc dy y yx dx y ②若Mx Ny M -??-??只含有y ，则原⽅程有积分因⼦.=-??-??dx Mxn y m e y )(µ，即0),()(),()(=+dy y x N y dx y x M y µµ是恰当的③若NxN y M ??-??只含y ，则?=??-??dy n xny m e y )(µ，即0),()(),()(=+dy y x N x dx y x M x µµ是恰当的④若MN xN y M -??-??,只含)(y x +，则?=++-??-??)()(y x d M N xny m e y x µ⑤若xMyN x N y M -??-??，只含有)(xy ，则?==??-)()(xy d xM yN x n y m e xy µ三、⼀阶微分⽅程的解的存在定理（1）研究的⽬的（2）解存在但不唯⼀的例⼦10,100)(22<≤<≤≤=-=?-=?=c x c x c x y c x y y dx dy其中（3）解的存在性定理⼀阶显⽰⽅程：),(y x f dxdy=……)1.3( 初值问题：==00)(),(y x y y x f dx dy ……)2.3(定理)1.1.3(存在唯⼀性定理如果)1.3(的),(y x f 在R ：b y y a x x ≤-≤-||,||00上满⾜：（1）在R 上连续（2）在R 上关于y 满⾜lipshit 条件，则初值问题)2.3(在区间h x x ≤-||0上上存在唯⼀解.其中),(y x f 对y 满⾜lipshit 条件是指，0>?L 常数,对R 中?两点),(),,.(1210y x y x 均有不等式成⽴：|||),(),(|2121y y L y x f y x f -≤-.20k y x y x f M mba h ∈=),(|),(|max ),,min( ⼏何解释：线段场定义)1.3(中的),(y x f 在2R k ∈内有定义，对R 中?点),(y x ，以),(y x 为中⼼，作⼀单位线段),(y x f k =，称为在点),(y x 的浅素。
2016春作业实验(1)常微分方程
1. 分别用Euler 法和ode45解下列常微分方程并与解析解比较：(1) ,(0)1, 03y x y y x '=+=<<function [ t,y ] = euler(f,ts,y0,h)t=ts(1):h:ts(2);y(1)=y0;for i=1:length(t)-1y(i+1)=y(i)+h*f(t(i),y(i));endt=t';y=y';endf=(t,y)t+y;[t1,y1]=euler(f,[0,3],1,0.05);[t2,y2]=ode45(f,[0,3],1);plot(t1,y1,'.-',t2,y2,'ro')hold ony3=dsolve('Dy=x+y','y(0)=1','x')ezplot(y3,[0,3])hold offlegend('euler','ode45','解析解');（2）22()5()3()45,(0)2,(0)1, 02tx t x t x t e x x t ''''--===<<f=(t,x)[2*x(2);5*x(2)+3*x(1)+45*exp(2*t)];[t1,y1]=ode45(f,[0,2],[2,1]);plot(t1,y1)2. 求一通过原点的曲线，它在(,)x y 处的切线斜率等于22,0 1.57.x y x +<<若x 上限增为1.58，1.60会发生什么？function dy = odefun_2(x,y)dy=2*x+y^2;dy=dy(:);end[t1,y]=ode45('odefun_2',[0,1.58],0) plot(t1,y);[t2,y]=ode45('odefun_2',[0,1.60],0) plot(t2,y);3. 求解刚性方程组：112121221000.25999.750.5,(0)1,050.999.751000.250.5,(0)1,y y y y x y y y y '=-++=⎧<<⎨'=-+=-⎩function Dy=fun(t,y)Dy=zeros(2,1);Dy(1)=-1000.25*y(1)+999.75*y(2)+0.5;Dy(2)=999.75*y(1)-1000.25*y(2)+0.5;[t,y]=ode15s('fun',[0,5],[1,-1]); plot(t,y(:,1),'o',t,y(:,2),'k-','LineWidth',2);4. (广告效应) 某公司生产一种耐用消费品，市场占有率为5%时开始做广告，一段时间的市场跟踪调查后，该公司发现：单位时间购买人口百分比的相对增长率与当时还没有买的百分比成正比，且估得此比例系数为0.5。
常微分方程数值解实验
多步法，Gear’s反向
数值积分，精度中等
若ode45失效时，
可尝试使用
ode23s
刚性
一步法，2阶Rosebrock算法，
低精度。
当精度较低时，
计算时间比ode15s短
odefx为显式常微分方程 中的 ，t为求解区间，要获得问题在其他指定点 上的解，则令t=[t0,t1,t2,…]（要求 单调），y0初始条件。
MATLAB 中有几个专门用于求解常微分方程的函数,它们的设计思想基于Runge-Kutta方法,基本设计思想为：从改进的欧拉方法比欧拉方法精度高的缘由着手,如果在区间[ x1, xi+1]多取几个点的斜率值,然后求取平均值,则可以构造出精度更高的计算方法。 这些函数主要包括：ode45、ode23、ode15s、ode113、ode23s、ode23t、ode23tb. 其中最常用的是函数ode45,该函数采用变步长四阶五阶Runge-Kutta法求数值解,并采用自适应变步长的求解方法。ode23采用二阶三阶Runge-Kutta法求数值解,与ode45类似,只是精度低一些。ode15s用来求刚性方程组。
43
4月22日
588
666
28
46
4月23日
693
782
35
55
4月24日
774
863
39
64
4月25日
877
954
42
73
4月26日
988
1093
48
76
4月27日
1114
1255
56
78
4月28日
1199
1275
59
78
4月29日
常微分方程实验报告
常微分方程实验报告《常微分方程》综合性实验实验报告实验班级05应数（3）学生姓名江晓荣学生学号200530770314指导老师方平华南农业大学理学院应用数学系实验微分方程在数学建模中的应用及数值解的求法一、实验目的1．了解常微分方程的基本概念。
2．常微分方程的解了解析解和数值解。
3．学习、掌握MA TLAB 软件有关求解常微分方程的解析解和数值解的有关命令。
4. 掌握微分方程在数学建模中的应用。
二、实验内容1．用MA TLAB 函数dsolve 符号求解常微分方程的通解和特解。
2．用MA TLAB 软件数值求解常微分方程。
三、实验准备1．用MA TLAB 求常微分方程的解析解的命令用MA TLAB 函数dsolve 求常微分方程()(,,,,,,)0n F x y y y y y ''''''= （7.1）的通解的主要调用格式如下：S=dsolve('eqn', 'var')其中输入的量eqn 是改用符号方程表示的常微分方程(,,,2,)0F x y Dy D y Dny = ，导数用D 表示，2阶导数用D2表示，以此类推。
var 表示自变量，默认的自变量为t 。
输出量S 是常微分方程的解析通解。
如果给定常微分方程（7.1）的初始条件()00010(),(),,()n n y x a y x a y x a '=== ，则求方程（7.1）的特解的主要调用格式如下：S=dsolve('eqn', 'condition1 ',…'conditonn ','var')其中输入量eqn ，var 的含义如上，condition1,…conditonn 是初始条件。
输出量S 是常微分方程的特解。
2．常微分方程的数值解法除常系数线性微分方程可用特征根法求解、少数特殊方程可用初等积分法求解外，大部分微分方程无解析解，应用中主要依靠数值解法。
电大《常微分方程》形成想考核作业参考答案.资料讲解
常微分方程第一、二、三次作业参考答案1、给定一阶微分方程2dyx dx=： （1） 求出它的通解；解：由原式变形得：2dy xdx =.两边同时积分得2y x C =+.（2） 求通过点（2，3）的特解；解：将点（2,3）代入题（1）所求的得通解可得：1C =-即通过点（2,3）的特解为：21y x =-.（3） 求出与直线23y x =+相切的解；解：依题意联立方程组：223y x Cy x ⎧=+⎨=+⎩故有：2230x x C --+=。
由相切的条件可知：0∆=，即2(2)4(3)0C --⨯-+=解得4C =故24y x =+为所求。
（4） 求出满足条件33ydx =⎰的解。
解：将2y x C =+代入330dy =⎰，可得2C =-故22y x =-为所求。
2、求下列方程的解。
１）3x y dydx-= 2）233331dy x y dx x y -+=-- 解：依题意联立方程组：23303310x y x y -+=⎧⎨-+=⎩解得：2x =，73y =。
则令2X x =-，73Y y =-。
故原式可变成：2333dY x y dX x y-=-. 令Yu X =，则dy Xdu udx =+，即有 233263u dxdu u u x-=-+. 两边同时积分，可得122(263)||u u C X --+= .将732y u x -=-，2X x =-代入上式可得： 12227()614323|2|2(2)y y C x x x -⎛⎫- ⎪--+=- ⎪-- ⎪⎝⎭.即上式为所求。
3、求解下列方程:1)24dyxy x dx+=. 解：由原式变形得：22dyxdx y=-. 两边同时积分得：12ln |2|y x C --=+. 即上式为原方程的解。
2)()x dyx y e dx-=. 解：先求其对应的齐次方程的通解：()0dyx y dx -=. 进一步变形得：1dy dx y=. 两边同时积分得：x y ce =.利用常数变异法，令()xy c x e =是原方程的通解。
数学实验——常微分方程数值解
实验4 常微分方程数值解分1 黄浩 43一、实验目的1.掌握用MATLAB软件求微分方程初值问题数值解的方法；2.通过实例学习用微分方程模型解决简化的实际问题；3.了解欧拉方法和龙格-库塔方法的基本思想和计算公式，及稳定性等概念。
二、实验内容1.《数学实验》第一版（问题2）问题叙述：小型火箭初始重量为1400kg，其中包括1080kg燃料。
火箭竖直向上发射时燃料燃烧率为18kg/s，由此产生32000N的推力，火箭引擎在燃烧用尽时关闭。
设火箭上升时空气阻力正比于速度的平方，比例系数为0.4kg/m，求引擎关闭瞬间火箭的高度、速度、加速度，及火箭到达最高点时的高度和加速度，并画出高度、速度、加速度随时间变化的图形。
模型转换及实验过程：（一）从发射到引擎关闭设火箭总质量为m，上升高度为h，瞬时速度为v，瞬时加速度为a，由燃料燃烧时间t=60s，可列如下的方程组：t∈[0,60]−9.8因此，上述方程为二元常微分方程组，选择t为自变量，h和v为因变量进行分析。
初值条件：h(0)=0 ,v(0)=0对上述模型，使用ode45()函数求数值解（程序见四.1、四.2），结果如下：由上表可知，引擎关闭瞬间，火箭的高度为12189.78m，速度为267.26m/s，加速度为0.9170m/s2，火箭至此已飞行60s而高度、速度、加速度随时间的变化曲线如下：(二)从引擎关闭到最高点设引擎关闭时，t1=0，由上一问的结果可知，h(0)=12189.78m，v(0)= 267.26m/s，m=320kg，则可列二元常微分方程组如下：9.8因此，可选择t1为自变量，h、v为因变量进行分析（程序见四.3、四.4），实验结果如下：由上表可知，当t1∈[11,12]时，v(t1)有零点，即该区间内某时刻火箭达到最高点。
再进行更细致的实验（程序略），设步长为0.01，观察该区间内v(t1)的零点，如下表所示：可以看出，当t1=11.30s，即总时间t=71.30s时，火箭达到最高点，高度为13115.36m，加速度为-9.8m/s2。
16秋华师《常微分方程》在线作业
B. y=x^2
C. y=e^(3x)
D. y=e^(2x)
正确答案：
9.按照微分方程通解定义，y''=sinx的通解是（）。
A. -sinx+C1x+C2
B. -sinx+C1+C2
C. sinx+C1x+C2
D. sinx+C1x+C2
正确答案：
10.方程组dY/dx=F(x,Y),x∈R,Y∈R^n的任何一个解的图象是（）维空间中的一条积分曲线.
B. dy/dx=k(x-a)(b-y)(k,a,b是常数)
C. dy/dx-siny=x
D. y'+xy=y^2*e^x
正确答案：
8.方程dy/dx=(1-y^2)(1/2)的常数解有（）。
A. y=1
B. y=-1
C. y=0
D. y为所有正数
正确答案：
9.方程y''+4y'+4y=0的基本解组有（）
C. s=-(1/2)g*t^2
D. s=1/2g*t^2
正确答案：
20.微分方程y'-y=0满足初始条件y(0)=1的特解为（）。
A. e^x
B. e^x-1
C. e^x+1
D. 2-e^x
正确答案：
华师《常微分方程》在线作业
二、多选题（共10道试题，共20分。）
1.下列哪些不可以作为变量可分离方程M(x)N(y)dx+p(x)q(y)dy=0的积分因子？（）
A. e^(2x)与2e^(2x)
B. e^(-2x)与xe^(-2x)
15春华师《常微分方程》在线作业答案
D. y'+y^2=0
正确答案：ABC
7.下列微分方程中，哪些不是可分离变量？（）。
A. dy/dx+y/x=e
B. dy/dx=k(x-a)(b-y)(k,a,b是常数)
C. dy/dx-siny=x
D. y'+xy=y^2*e^x
正确答案：ACD
8.在下列函数中，不可能是微分方程y''+y=0的解的函数有（）。
D. y''+y=cosx
正确答案：C
19. f(y)连续可微是保证方程dy/dx=f(y)解存在且唯一的（）条件．
A.必要
B.充分
C.充分必要
D.必要非充分
正确答案：B
20.方程dy/dx=y^(1/2)+1（）奇解．
A.有一个
B.有两个
C.无
D.有无数个
正确答案：C
华师《常微分方程》在线作业
二、多选题（共10道试题，共20分。）
1.下列哪些不可以作为变量可分离方程M(x)N(y)dx+p(x)q(y)dy=0的积分因子？（）
A. 1/(N(y)+P(x))
B. 1/(N(y)-P(x))
C. 1/(N(y)P(x))
D. 1/(P(x)-N(y))
正确答案：ABD
2.下列函数中，哪些不是微分方程xdy+ydy=0的通解？（）。
C. e^x+1
D. 2-e^x
正确答案：A
13.微分方程y''-4y'+4y=0的两个线性无关解是（）。
A. e^(2x)与2e^(2x)
常微分方程课后作业
习题2.11.xy dxdy 2=,并求满足初始条件：x=0,y=1的特解.解：对原式进行变量分离得。
故它的特解为代入得把即两边同时积分得：e e xx y c y x x c y c y xdx dy y22,11,0,ln ,212=====+==17.yy y x x xy x dxdy -+++=3232332解：原方程化为123132;;;;;)123()132(2222222222-+++=-+++=y x y x dxdy y x y y x x dxdy令)1.......(123132;;;;;;;;;;;;,22-+++===u v u v dvdu v x u y 则方程组，，，）；令，的解为（111101230132+=-=-⎩⎨⎧=-+=++u Y v Z u v u v则有⎪⎪⎩⎪⎪⎨⎧++==+=+z y zy dzdy y z y z 23321023032）化为，，，，从而方程（，令 )2.( (232223322)，，，，，所以，，则有ttdzdt ztt dzdt zt dzdt zt dzdy zy t +-=++=++==当是原方程的解或的解。
得，是方程时，，即222222)2(1022x yx yt t-=-=±==-当c x y xy dz zdt tt t5222222)2(12223022+-=+=-+≠-两边积分的时，，分离变量得另外cx y xy x yx y522222222)2(2+-=+-=-=原方程的解为，包含在其通解中，故，或4．dxdy nxx e y nx =-， n 为常数.解：原方程可化为：dxdy nx x e y nx +=)(c dx ex e ey dxx nnx dxx n+⎰⎰=⎰-)(c e x xn+= 是原方程的解. 6．dx dy 234xyx x +=解：dxdy 234xyx x +==23yx +xy令xy u = 则 ux y =dxdy =u dxdu x+因此：dxdu xu +=2ux ，21udxdu=，dx du u =2，c x u +=331c x x u +=-33 （*） 将xy u =带入 （*）中 得：3433cx x y =-是原方程的解.15331dy dxxy x y=+33dx yx y x dy=+这是n=3时的伯努利方程。
数学实验第二次作业——常微分方程数值求解
实验4常微分方程数值解实验目的：1.练习数值积分的计算；2.掌握用MATLAB软件求微分方程初值问题数值解的方法；3.通过实例学习用微分方程模型解决简化的实际问题；4.了解欧拉方法和龙格——库塔方法的基本思想和计算公式，及稳定性等概念。
实验内容：3.小型火箭初始质量为1400kg，其中包括1080kg燃料，火箭竖直向上发射是燃料燃烧率为18kg/s，由此产生32000N的推力，火箭引擎在燃料用尽时关闭。
设火箭上升是空气阻力正比于速度的平方，比例系数为0.4kg/m，求引擎关闭瞬间火箭的高度，速度，加速度，及火箭到达最高点是的高度，速度和加速度，并画出高度，速度，加速度随时间变化的图形。
解答如下：这是一个典型的牛顿第二定律问题，分析火箭受力情况；先规定向上受力为正数建立数学模型：A燃料未燃尽前，在任意时刻（t<60s）火箭受到向上的-F=32000N，向下的重力G=mg，g=9.8，向下的阻力f=kv^2, k=0.4, v表示此时火箭速度；此时火箭收到的合力为F1=（F-mg-f）；火箭的初始质量为1400kg，燃料燃烧率为-18kg/s;此刻火箭质量为m=1400-18*t根据牛顿第二定律知，加速度a=F1/m=(F-mg-f)/(m-r*t)=(32000-(0.4.*v.^2)-9.8.*(1400-18.*t))由此可利用龙格-库塔方法来实现，程序实现如下Function [dx]=rocket[t,x] %建立名为rocket的方程m=1400;k=0.4;r=-18;g=9.8; %给出题目提供的常数值dx=[x(2);(32000-(k*x(2)^2)-g*(m+r*t))/(m+r*t)];%以向量的形式建立方程[a]=(32000-(k*x(2)^2)-g*(m+r*t))/(m+r*t); %给出a的表达式End;ts=0:60; %根据题目给定燃烧率计算出燃料燃尽的时间，确定终点x0=[0,0]; %输入x的初始值[t,x]=ode15s(@rocket,ts,x0); %调用ode15s计算[t,x];h=x(:,1);v=x(:,2);plot(t,x(:,1)),grid; %绘出火箭高度与时间的关系曲线title('h-t');xlabel('t/s');ylabel('h/m'),pause;plot(t,x(:,2)),grid ; %绘出火箭速度与时间的曲线关系title('v-t');xlabel('t/s');ylabel('v/m/s'),pause;a=(32000-(0.4.*v.^2)-9.8.*(1400-18.*t))/(1400-18.*t); plot(t,a),grid; %绘出火箭加速度与时间的曲线关系title('a-t');xlabel('t/s'),ylabel('a/m^2/s'),pause火箭高度随时间变化的曲线火箭速度随时间变化的曲线火箭加速度随时间变化的曲线数据过多，故截取部分如下第一列为时间，第二列为火箭高度，第三列为火箭速度由此可以，在t=60s时，即火箭燃料燃尽瞬间，引擎关闭瞬间，火箭将到达12912m的高度，速度为267,29m，加速度a=0.9m/s^2B燃料燃尽之后，与A 类似，分析受力如下火箭受到向上的F=0向下的重力G=mg，g=9.8，向下的阻力f=kv^2, k=0.4, v表示此时火箭速度；此时火箭收到的合力为F2=（-mg-f）；火箭的初始质量为320kg，恒定根据牛顿第二定律,加速度a=F2/m=-g-0.4v^2/320;程序实现如下function [ dx ] = rocket2( t,x ) %建立以rocket2为名的函数dx=[x(2);-9.8-0.4.*x(2).^2/320]; %以向量的形式建立方程ts=60:120; %给出初始时刻，估计终点时刻x0=[12190,267.26]; %给出x初始值[t,x]=ode15s(@rocket2,ts,x0); %调用ode15s计算[t,x]plot(t,x(:,1)),grid; %绘出火箭高度随时间变化的曲线title('h-t');xlabel('t/s'),ylabel('h/m'),pause;plot(t,x(:,2)),grid; %绘出火箭速度随时间的变化曲线title('v-t');xlabel('t/s'),ylabel('v/m/s'),pause;v=x(:,2);a=-9.8-0.4*v.^2/320; %给出加速度的具体表达式plot(t,a),grid; %绘出火箭加速度随时间变化的曲线title('a-t');xlabel('t/s'),ylabel('a/m^2/s'),pause得到的曲线图形如下火箭高度随时间的变化曲线从图中可以大致看出，最高点在13km左右，火箭速度随时间的变化曲线加速度随时间变化曲线如下数据表格大致如下从图表中可以看出，在71s左右速度到达0，即此时到达最高处，高度为13117m加速度为-9.8m/m/s^2;本题总结：这道题是典型的物理牛顿力学的题目，通过受力的正确分析，可以知道，以[h,v]为向量建立微分方程即可求解，h的微分是速度v，速度v的微分是加速度a解题过程中存在的难点是：取值步长不太容易确定，而且是哪种算法不确定，先用ode15s 速度较快，ode23s速度差不太多，其他两种速度较慢，等待时间较长5．一只小船渡过宽为d 的河流，目标是起点A 正对着的另一岸B 点。
奥鹏东师 《常微分方程》练习题答案.docx
《常微分方程》练习题一参考答案练习题第1套参考答案 一. 填空题1、全平面.2、1,1x y =-=-3、3y Cx C =+ 4、线性无关,(或朗斯基行列式不等于零) 5、开二. 单项选择题1.A,2.C,3.B,4.C,5.B三. 简答题1.0y >时对应通解是2(),.4x C y C x +=-≤<∞ 0y <时对应通解是2(),.4x C y x C +=--∞≤<- 2.是.四. 计算题 1、通积分为1x y Ce y -=. 2、通解为411().4y C x x =+ 3、通积分为21.x y C y += 4、通解为121cos sin cos .2x C t C t t t =+- 5、通解为27124151t t x C e C e y -⎡⎤⎡⎤⎡⎤=+⎢⎥⎢⎥⎢⎥-⎣⎦⎣⎦⎣⎦五. 应用题1. 设物体在t 时刻的下落速度为().v v t =在t 时刻物体所受的力,f mg kv =-k 为阻力系数,由牛顿第二运动定律,得方程dv m mg kv dt =- 即 ()dv k mg v dt m k=-- 解得 kt mmg v Ce k -=+ 代入初值条件(0)0v =, 得初值解 ()(1)kt m mg v t e k -=- 令t →+∞,得极限速度1.mg v k=2. 证明:因为0x 在取极值有1020()()0y x y x ''== 此时12(),()y x y x 的朗斯基行列式在0x 点的值为 1020102001020()()()()()0()()0y x y x y x y x W x y x y x ==='' 所以, 12(),()y x y x 不能为基本解组.练习题第2套参考答案 一、填空题1、(,)-∞+∞.2、0y >的右半平面3、,0,1,2,y k k π==±±L4、 22,xx exe -- 5、n二、单项选择题1.B,2.A,3.D,4.C,5.D三、简答题化成等价积分方程，用逐次逼近法求积分方程解。
常微分方程-习题作业-第六章第四节作业及详细解答
dv dt
=
√ − 2u − v,
由此不难画出变换后的系统在 uv 平面上的相图. 由此可画出原系统在 xy 平面上的相图. (4) 容易求得平衡点为 (−1, −1). 引进平移变换 u = x + 1, v = y + 1 可将系统化为
du dt
=
5u + 3v,
dv dt
=
−3u − 5v.
其系数矩阵的特征值为λ1 = 4, λ2 = −4, 是一对相异实根, 符号相反, 因此平衡点 (−1, −1) 为原系统的鞍点, 不稳定. 它有两个特殊方向, 容易求得对应于 λ1 的特征向量 ξ1 = (3, −1)T , 对应于 λ2 的特征向量 ξ2 = (1, −3)T , 相应于 ξ1 的直线上面的轨道都是继续沿着 它背离平衡点, 相应于 ξ2 的直线上面的轨道都是继续沿着它趋向平衡点, 由此可画出原系 统在 xy 平面上的相图.
5. 引入极坐标并画出下面系统的相图：
dx dt
=
x(x2
+
y2
−
1),
dy dt
=
y(x2
+ y2 − 1).
解: 令 x = r cos θ, y = r sin θ, 原系统化成:
dr dt=r来自r2− 1),dθ dt
=
0.
易知它有特解 (r(t), θ(t)) ≡ (0, t0) 及 (r(t), θ(t)) ≡ (1, t0), 其中 t0 为任意常数. 它们都对应 着原系统的奇点. 因此原系统有奇点 (0, 0) 及奇线 x2 + y2 = 1. 容易由极坐标方程看出除奇 点及奇线外该方程组的轨线族为相平面上的一族射线 θ(t) = t0, 在奇线 x2 + y2 = 1 内部, 它趋于奇点, 在奇线 x2 + y2 = 1 外部, 它远离奇线 x2 + y2 = 1. 故原点是稳定的星形结点.
微分方程数值方法大作业一
微分方程数值方法大作业一、操作任务给定初值问题22x u t u ∂∂=∂∂ 0<x<1,0<t<0.1u(0,t)=u(1,t)=0 0<t<0.1u(x,0)=sinπx 0≤x≤1研究上述问题Grank-Nicholson 格式的稳定性二、差分格式先空间方向离散：取x=x j 得： ,即得半离散格式：记u(t)=[u(x 1,t), u(x 2,t) , …… u(x N-1,t)]TF(t) =[f(x 1,t)], f(x2,t)] ……f(x N-1,t)]TN h 1=h j x j ⨯=，，j=0,1,2……N),(22),(t x x u t x t u j j ∂∂=∂∂[]),(),(2),(1),(112t x u t x u t x u h t x t u j j j j -++-=∂∂j=0,1,2……N21112112---=A格式变为：U(0)=[)()(),(121-N x x x ϕϕϕ ]T用梯形格式——Grank-Nicholason 格式则格式变为：)0(0u u =三、程序function x=EqtsForwardAndBackward(L,D,U,b)%追赶法求解三对角线性方程组Ax=b%x=EqtsForwardAndBackward(L,D,U,b)%x:三对角线性方程组的解%L:三对角矩阵的下对角线，行向量%D:三对角矩阵的对角线，行向量%U:三对角矩阵的上对角线，行向量%b:线性方程组Ax=b 中的b ，列向量n=length(D);m=length(b);n1=length(L);n2=length(U);if n-n1 ~= 1 || n-n2 ~= 1 || n ~= m %检查参数的输入是否正确disp('输入参数有误！')x=' ';return ;end%追的过程 )(1)(2t Au h t dt du =⎥⎦⎤⎢⎣⎡+=-++)()(211)()(121k k k k t u t u A h t u t u τkk u rA I u rA I )21()21(1+=-+for i=2:nL(i-1)=L(i-1)/D(i-1);D(i)=D(i)-L(i-1)*U(i-1);endx=zeros(n,1);x(1)=b(1);for i=2:nx(i)=b(i)-L(i-1)*x(i-1);end%赶的过程x(n)=x(n)/D(n);for i=n-1:-1:1x(i)=(x(i)-U(i)*x(i+1))/D(i);endreturn;function [U x t] = PDEParabolicCN(uX,uT,M,N) %Crank-Nicolson隐式格式求解抛物型偏微分方程%输入参数：uX -空间变量x的取值上限% ? ?? ?? uT -时间变量t的取值上限% ? ?? ?? M -沿x轴的等分区间数% ? ?? ?? N -沿t轴的等分区间数h=uX/M;%计算空间方向步长tao=uT/N;%计算时间方向步长x=(0:M)*h;t=(0:N)*tao;r=h/(tao*tao);%网格比Diag=zeros(1,M-1);%矩阵的对角线元素Low=zeros(1,M-2);%矩阵的下对角线元素Up=zeros(1,M-2);%矩阵的上对角线元素for i=1:M-2Diag(i)=1+r;Low(i)=-r/2;Up(i)=-r/2;endDiag(M-1)=1+r;%计算初值和边值U=zeros(M+1,N+1);for i=1:M+1U(i,1)=sin(pi*x(i));endfor j=1:N+1U(1,j)=0;U(M+1,j)=0;endB=zeros(M-1,M-1);for i=1:M-2B(i,i)=1-r;B(i,i+1)=r/2;B(i+1,i)=r/2;endB(M-1,M-1)=1-r;%逐层求解，需要使用追赶法（调用函数EqtsForwardAndBackward）for j=1:Nb1=zeros(M-1,1);b1(1)=r*(U(1,j+1)+ U(1,j))/2;b1(M-1)=r*(U(M+1,j+1)+U(M+1,j))/2;b=B*U(2:M,j)+b1;U(2:M,j+1)=EqtsForwardAndBackward(Low,Diag,Up,b); endU=U';disp('请输入uX')uX=input('uX=');disp('请输入uT')uT=input('uT=');disp('请输入M')M=input('M=');disp('请输入N')N=input('N=');[U x t] = PDEParabolicCN(uX,uT,M,N);mesh(x,t,U);title('Grank-Nicholson，隐式格式，一维热传导方程解的图像') xlabel('空间变量 x')ylabel('时间变量t')zlabel('一维热传导方程的解 x')return;四、数据结果M=50时：M=100时：五、结论该方程组Grank-Nicholason格式稳定。
《常微分方程》练习题库参考答案
华中师范大学网络教育学院 《常微分方程》练习测试题库参考答案一、判断说明题1、在线性齐次方程通解公式中C 是任意常数而在常数变易法中C （x ）是x 的可微函数。
将任意常数C 变成可微函数C （x ），期望它解决线性非齐次方程求解问题，这一方法成功了，称为常数变易法。
2、因p(x)连续，y(x)= y 0exp(-dx x⎰0x p(x))在p(x)连续的区间有意义，而exp(-dx x⎰x p(x))＞0。
如果y 0＝0，推出y(x)=0,如果y(x)≠0,故零解y(x)=0唯一。
3、（1） 它是常微分方程，因为含有未知函数的导数，f,g 为已知函数，y 为一元函数，所建立的等式是已知关系式。
（2） 它是常微分方程，理由同上。
（3） 它不是常 微分方程，因y 是未知函数，y(y(y(x)))也是未知的，所建立的等式不是已知关系式。
4、微分方程求解时，都与一定的积分运算相联系。
因此，把求解一个微分方程的过程称为一个微分方程。
微分方程的解又称为（一个）积分。
5、 把微分方程的通解用初等函数或通过它们的积分来表达的方法。
注意如果通解能归结为初等函数的积分表达，但这个积分如果不能用初等函数表示出来，我们也认为求解了这个微分方程，因为这个式子里没有未知函数的导数或微分。
6、 y `＝f(x,y)主要特征是f(x,y)能分解为两个因式的乘积，其中一个因式仅含有x,另一因式仅含y ，而方程p(x,y)dx+q(x,y)dy=0是可分离变量方程的主要特征，就像f(x,y)一样，p,q 分别都能分解成两个因式和乘积。
7、二元函数f(x,y)满足f(rx,ry)=r mf(x,y),r.>0,则称f(x,y)为m 次齐次函数。
m=0则称它为0次齐次函数。
8、如果f(x,y)是0次齐次函数，则y `＝f(x,y)称为齐次方程。
如果p(x,y)和q(x,y)同为m 次齐次函数，则pdx+qdy=0为齐次方程。
《常微分方程》网络作业
《常微分方程》网络作业1作业1给定一阶微分方程2dyx dx=： （1）求出它的通解；（2）求通过点（2，3）的特解； （3）求出与直线23y x =+相切的解； （4）求出满足条件33ydx =⎰的解。
解：（1）2y x c =+（其中c 为任意常数） （2）把点（2，3）代入2y x c =+得c ＝﹣1 ∴ 所求通过点（2，3）的特解为21y x =- （3）∵ 所求的解与直线23y x =+相切 ∴22dyx dx==，得x ＝1 把x ＝1代入23y x =+得y ＝5这表明所求的解与直线23y x =+相切于点（1，5） ∴ c ＝4∴ 所求的解为24y x =+ （4）把2y x c =+代入33ydx =⎰得320()3x c dx +=⎰，即331()33x cx +=，即1(273)033c ⨯+-=，解得c ＝﹣2∴ 所求的解为22y x =- 作业2求下列方程的解： 1）3x y dydx-= 2）233331dy x y dx x y -+=--解：1）3x y dydx-=这是一个变量分离方程 变量分离得33yxdy dx =两边同时积分得33yxc =+（其中c 为任意常数）2）由23303310x y x y -+=⎧⎨--=⎩解得4113x y =⎧⎪⎨=⎪⎩，作变换4113x y ξη=-⎧⎪⎨=-⎪⎩，则有2333d d ηξηξξη-=- 令u ηξ=，则有226333du u u d u ξξ-+=-，变量分离，并两边同时积分得22(263)u u c ξ-+= 把u ηξ=及4113x y ξη=-⎧⎪⎨=-⎪⎩代入可得2211112(4)6(4)()3()33x x y y c ----+-=即22126362x xy y x y c -+++=（其中c ，c 1为任意常数），这即为所求的通解。
作业3 求解下列方程:1)24dyxy x dx += 2)()x dyx y e dx -=53)dy y xy dx-=2234)42(1)0x y dx x y dy +-=解：1）该方程可化为2(2)dyx y dx=-，这是变量分离方程 变量分离得22dyxdx y=- 两边积分得原方程的通解为22x y ce -=+（c 为任意常数）2）齐线性方程dyy dx=的通解为x y ce =（c 为任意常数） 设()xy C x e =为方程的解，则()1dC x dx x=，解得()ln C x x c =+（c 为任意常数） ∴ 原方程的通解为(ln )xy e x c =+（c 为任意常数） 3）这是n ＝5的贝努利方程。
计算方法--常微分方程求解实验
实验五 常微分方程求解实验一、 实验目的通过本实验学会对给定初值我呢他，用欧拉法、改进欧拉法、四阶龙格-库塔法求数值解和误差，并比较优缺点.对给定刚性微分方程，求其数值解，并与精确解比较，分析计算结果.二、 实验题目1. 解初值问题各种方法比较 实验题目：给定初值问题⎪⎩⎪⎨⎧=≤<+=,0)1(,21,e d d y x x xyx y x 精确解为)e -e (xx y =，按 （1） 欧拉法，步长;1.0,025.0==h h （2） 改进欧拉法，步长;01.0,05.0==h h （3） 四阶标准龙格-库塔法，步长;1.0=h求在节点)10,,2,1(1.01 =+=k k x k 处的数值解及误差，比较各方法的优缺点. 2. 刚性方程计算实验题目：给定刚性微分方程⎪⎩⎪⎨⎧=≤<-+-=-,2)0(,50,600e 8.1199600d d 1.0y x y xyx 其精确解为12e e )(-0.1x -600x-+=x y .任选取一显示方法，取不同的步长求解，并分析计算结果.三、 实验原理1.欧拉格式由数值微分的向前差商公式可以解决初值问题（6.1）()()⎩⎨⎧=≤≤=00',,,y x y b x a y x f y 中的导数的数值计算问题：()()()()(),'1h x y x y h x y h x y x y n n n n -=-+≈+由此可得()()().'1n n n x hy x y x y +≈+（6.1）实际上给出()()()()()().,','n n n x y x f x y x y x f x y =⇒=于是有()()()().,1n n n n x y x hf x y x y +≈+再由()()11,++≈≈n n n n x y y x y y 得().,1,0,,111 =+=+++n y x hf y y n n n n （6.2） 递推公式（6.2）称为欧拉格式。
数值分析-常微分数值解的求法C语言
本科生课程设计报告实习课程数值分析学院名称管理科学学院专业名称信息与计算科学学生姓名学生学号指导教师实验地点实验成绩二〇一六年六月二〇一六年六摘要常微分方程数值解法是计算数学的一个分支.是解常微分方程各类定解问题的数值方法.现有的解析方法只能用于求解一些特殊类型的定解问题，实用上许多很有价值的常微分方程的解不能用初等函数来表示，常常需要求其数值解.所谓数值解，是指在求解区间内一系列离散点处给出真解的近似值.这就促成了数值方法的产生与发展.关键词：数值解法；常微分1. 实验内容和要求常微分方程初值问题cos22sin 22 0<<2(0)1xy y x x xe x y -'⎧=-+-+⎨=⎩有精确解2()cos(2)x y x x e x -=+。
要求：分别取步长h = 0.1，0.01，0.001，采用改进的Euler 方法、4阶经典龙格-库塔R －K 方法和4阶Adams 预测-校正方法计算初值问题。
以表格形式列出10个等距节点上的计算值和精确值，并比较他们的计算精确度。
其中多步法需要的初值由经典R-K 法提供。
运算时，取足以表示计算精度的有效位。
2. 算法说明2.1函数及变量说明表1 函数及变量定义1、欧拉方法：1()()(,())i i i i y x y x hf x y x +=+ （i=0，1，2，3，......n-1） (0)y a= (其中a 为初值)2、改进欧拉方法：~1~111()()(,())()()[(,())(,())]2(0)i i i i i i i i i i y x y x hf x y x hy x y x f x y x f x y x y a ++++=+=++=（i=0，1，2......n-1） (其中a 为初值)3、经典K-R 方法： 11213243()6(,)(,)22(,)22(,)i i i i i i i i i i h y y K f x y h h K f x y K h h K f x y K K f x h y hK +⎧=+⎪⎪=⎪⎪⎪=++⎨⎪⎪=++⎪⎪=++⎪⎩4、4阶adams 预测-校正方法预测：校正：Adsms 内插外插公式联合使用称为Adams 预测-校正系统，利用外插公式计算预测，用内插公式进行校正。
2013春数学实验基础 实验报告(1)常微分方程
实验一 常微分方程1. 分别用Euler 法和ode45解下列常微分方程并与解析解比较： (1)30,1)0(,<<=+='x y y x y编写Euler 法的matlab 函数，命名为euler.mfunction [t,y]=euler(odefun,tspan,y0,h) t=tspan(1):h:tspan(2); y(1)=y0;for i=1:length(t)-1y(i+1)=y(i)+h*feval(odefun,t(i),y(i)); end t=t'; y=y';下面比较三者的差别： % ode45odefun=inline('x+y','x','y'); [x1,y1]=ode45(odefun,[0,3],1); plot(x1,y1,'ko');pause hold on ; % Euler·¨[x2,y2]=euler(odefun,[0,3],1,0.05); plot(x2,y2,'r+');pause hold on ;% 解析解y0=dsolve('Dy=t+y','y(0)=1'); ezplot(y0,[0,3]);pause hold off ;legend('ode45','euler 法','解析解');实验一 常微分方程Euler 法只有一阶精度，所以实际应用效率比较差，而ode45的效果比较好，很接近真实值。
(2) 20.01()2sin(),(0)0,(0)1,05y y y t y y t ''''-+===<<先写M 文件ex1_2fun.mfunction f=ex1_2fun(t,y) f(1)=y(2);f(2)=0.01*y(2).^2-2*y(1)+sin(t); f=f(:); % ode45[t1,y1]=ode45(@ex1_2fun,[0,5],[0;1]); plot(t1,y1(:,1),'ko');% 解析解s=dsolve('D2y-0.01*(Dy)^2+2*y=sin(t)','y(0)=0','Dy(0)=1','t')s =[ empty sym ]%由此可知该微分方程无解析解2. 求一通过原点的曲线，它在(,)x y 处的切线斜率等于22,0 1.57.x y x +<<若x 上限增为1.58，1.60会发生什么？odefun=inline('2*x+y^2','x','y'); subplot(1,4,1);[x1,y1]=ode45(odefun,[0,1.57],0); plot(x1,y1,'r*'); title('上限1.57'); subplot(1,4,2);[x2,y2]=ode45(odefun,[0,1.58],0); plot(x2,y2,'bo'); title('上限1.58');实验一 常微分方程subplot(1,4,3);[x3,y3]=ode45(odefun,[0,1.6],0); plot(x3,y3,'k'); title('上限1.60'); subplot(1,4,4);plot(x1,y1,'r*');hold on ; plot(x2,y2,'bo');hold on ; plot(x3,y3,'k');hold off ;legend('上限1.57','上限1.58','上限1.60');上限1.5713上限1.5814上限1.6014结论：随着x 上界的增加，解趋于无穷大。
《常微分方程》作业参考答案
《常微分方程》作业参考答案一．求解下列方程1．x c y cos =2．通解为：x x c y sin cos +=3．dx x x dy 122-= ⎰⎰--=122)1(xx d dy 2ln 1y x c =-+ 1)0(==c y 2ln |1|1y x ∴=-+4．'(1)ln(1)y yyy x x x -=++ 令 xuy x yu =⇔= (1)ln(1)dyduu x u u u dx dx ∴=+=+++故 (1)ln(1)dux u u dx =++(1)ln(1)du dx u u x =++ ln(1)ln(1)d u dxu x +=+ln ln(1)ln ln u x c ∴+=+ ln(1)u cx +=cx e u =+1 cx e x y=+∴1 )1(-=cx e x y5. 可分离变量方程，通解为.)1)(1(222cx y x =++6．齐次方程，通解为 c x x yx y =++ln 422sin .7．全微分方程，通解为 .64224c y y x x =+-8．.0222=++y dx dyx dx y d9. 解为 .)3(3x x y -=10. 通解为 .2sin 222c y x y x =++11．方程为 .011222=+-y x dx dyx dx y d12．通解为 ).tan(21c x c y +=二．1．通解为：c e e x y +=2212. 通解为： t t e c c e c z y x 2321123101210⎪⎪⎪⎭⎫ ⎝⎛-+⎪⎪⎪⎭⎫ ⎝⎛+⎪⎪⎪⎭⎫ ⎝⎛-=⎪⎪⎪⎭⎫ ⎝⎛-3.0)0(0==y y 2121x y =52220121x x y += 4. x uN y uM ∂∂=∂∂ xu N x N u y u M y M u ∂∂+∂∂=∂∂+∂∂ 令 u y x =+22 y u d u d y u 2⋅=∂∂∴ x ud u d x u 2⋅=∂∂ u d u d x x N u u d u d y y M u 22+∂∂=+∂∂ ⎪⎪⎭⎫ ⎝⎛∂∂-∂∂=-∴y M x N u u d u d x y )(2故满定充要条件的表达式为：)(22y x xy y M xN +=--∂∂∂∂ϕ 5.)(2122y x v +=)(*dtdv)(22s x +-≤∠0 022≠+s x ∴（0.0）渐近稳定 6.一次近似方程为：⎪⎩⎪⎨⎧+=--=y x dtdy y x dt dx 32 特征方程为：012=++λλ 3-=∴∆＜0 P ＝1＞0 ∴0)Re(0)Re(21<<λλ, 则（0.0）局部渐过稳定. 7.01032=--λλ 5,221=-=λλx B x B x A x A y o 2sin )(2cos )(101*1+++=为x x y y y 2cos 10'3"=-- 之特解,±2λ不是特征根5=a 是特征方程的单根 x o e c x c x c x y 52122)(++=∴*故其通解为： 215221y y e c ec y x x +++=-8.特征根为：2.1.1321==-=λλλ 11-=λ所属的特征向量为：⎪⎪⎪⎭⎫ ⎝⎛-=532α12=λ所属的特征向量为：⎪⎪⎪⎭⎫ ⎝⎛=111β13=λ所属的特征向量为：γ⎪⎪⎪⎭⎫ ⎝⎛=101通解为：t t t e c e c e c z y x 2321101111531⎪⎪⎪⎭⎫ ⎝⎛+⎪⎪⎪⎭⎫ ⎝⎛+⎪⎪⎪⎭⎫ ⎝⎛-=⎪⎪⎪⎭⎫ ⎝⎛-9.0:)0(=o y y 2121x y =52220121x x y -= 10.特征方程为：01072=++λλ07>=p 010>=g 0>∆故 (0.0)为稳定结点11.1.一次近似方程为：⎪⎩⎪⎨⎧-=--=yx y x t d y d dt x d 0222=++∴λλ0)Re(1<λ 0)Re(2<λ ∴（0.0）为局部渐近稳定 2.)(2122y x v +=. )1)((2222)(-++=*y x y x l dt dv 故122<+y x 0<∴dtdv 故（0.0）局部渐近稳定. 12. 1.,00=y ,31),(3020001x dx x dx y x f y y x x==+=⎰⎰ .63131)91(),(730620102x x dx x x dx y x f y y x x+=+=+=⎰⎰ 2. ,),(22y x y x f += ∴ ,5),(max ),(==∈y x f M Dy x ,42max max ),(),(L y y f D y x D y x ===∂∂∈∈ .5252,1min ,min =⎭⎬⎫⎩⎨⎧=⎭⎬⎫⎩⎨⎧=m b a h则 .7564)52(32145)()(322=⋅⋅⋅≤-x y x y 13. 系数阵为 ,110111110⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡- 特征方程为 .0)1()det(2=--=-λλλE A E A λ-的初等因子为 2)1(,-λλ，通解为.101010101112321t t e t c e c c z y x ⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡⎪⎪⎪⎭⎫ ⎝⎛+⎪⎪⎪⎭⎫ ⎝⎛+⎪⎪⎪⎭⎫ ⎝⎛+⎪⎪⎪⎭⎫ ⎝⎛-=⎪⎪⎪⎭⎫ ⎝⎛14.证：设 [).),0()(..,0+∞∈∀≤>∃x M x f t s M .则[)+∞∈∀,0x ，有 .)1()(0)(0000M y e M y ds e Me y x y x x xx s x+≤-+=+≤--⎰[]),,0()(0x C x y ∈ ∴ [].,0,)(..,00x x M x y t s M ∈≤>∃令 {},,max 0M y M K += ∴ [).,0,)(+∞∈∀≤x K x y15.通解为 .)21(221xx e x x x c e c y -++=16．,2=α 特解为 ,1x y = 通解为 ).ln 21(221x x x c x c y +-+=。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
1、 分别用Euler 法与ode45解下列常微分方程并与解析解比较:
(1) ,(0)1, 03y x y y x '=+=<<
function [ t,y ] = euler(f,ts,y0,h)
t=ts(1):h:ts(2); y(1)=y0;
for i=1:length(t)-1
y(i+1)=y(i)+h*f(t(i),y(i)); end t=t'; y=y'; end
f=@(t,y)t+y;
[t1,y1]=euler(f,[0,3],1,0、05); [t2,y2]=ode45(f,[0,3],1); plot(t1,y1,'、-',t2,y2,'ro') hold on
y3=dsolve('Dy=x+y','y(0)=1','x') ezplot(y3,[0,3]) hold off
legend('euler','ode45','解析解');
（2）22()5()3()45,(0)2,(0)1, 02t
x t x t x t e x x t ''''--===<<
f=@(t,x)[2*x(2);5*x(2)+3*x(1)+45*exp(2*t)]; [t1,y1]=ode45(f,[0,2],[2,1]); plot(t1,y1)
2. 求一通过原点的曲线,它在(,)x y 处的切线斜率等于2
2,0 1.57.x y x +<<若x 上限增为1、58,1、60会发生什么？
function dy = odefun_2(x,y) dy=2*x+y^2; dy=dy(:); end
[t1,y]=ode45('odefun_2',[0,1、58],0) plot(t1,y);
[t2,y]=ode45('odefun_2',[0,1、60],0) plot(t2,y);
3、 求解刚性方程组:
11212
1221000.25999.750.5,(0)1,050.
999.751000.250.5,(0)1,y y y y x y y y y '=-++=⎧<<⎨
'=-+=-⎩
function Dy=fun(t,y) Dy=zeros(2,1);
Dy(1)=-1000、25*y(1)+999、75*y(2)+0、5; Dy(2)=999、75*y(1)-1000、25*y(2)+0、5; [t,y]=ode15s('fun',[0,5],[1,-1]);
plot(t,y(:,1),'o',t,y(:,2),'k-','LineWidth',2);
4、(广告效应) 某公司生产一种耐用消费品,市场占有率为5%时开始做广告,一段时间的市场跟踪调查后,该公司发现:单位时间内购买人口百分比的相对增长率与当时还没有买的百分比成正比,且估得此比例系数为0、5。
(1) 建立该问题的数学模型,并将解析解与数值解,并作以比较;
y’=0、5(1-y)
y=desolve('Dy=0、5-0、5*y','y(0)=0、05')
odefun=@(t,y)0、5-0、5*y;
[t1,y1]=ode45(odefun,[0,10],0、05);
t2=0:0、1:10;
y2=1-(19*exp(-t2/2))/20;
plot(t1,y1,'o',t2,y2,'k');
(2)厂家问:要做多少时间广告,可使市场购买率达到80%？
1-(19*exp(-t/2))/20=0、8
5、(肿瘤生长) 肿瘤大小V生长的速率与V的a次方成正比,其中a为形状参数,0≤a≤1;而其比例系数K随时间减小,减小速率又与当时的K值成正比,比例系数为环境参数b。
设某肿瘤参数a=1, b=0、1, K的初始值为2,V的初始值为1。
问
(1)此肿瘤生长不会超过多大？
k’=-bk,v’=k*v^a,得k’=-0、1k,v’=kv,且k(0)=2,v(0)=1,
[k,v]=dsolve('Dk=-0、1*k','Dv=k*v','k(0)=2','v(0)=1','t');
t=0:0、1:100;
v=exp(20)*exp(-20*exp(-t/10));
plot(t,v);
(2)过多长时间肿瘤大小翻一倍？
exp(20)*exp(-20*exp(-t/10))=2
(3)何时肿瘤生长速率由递增转为递减？
v’与v的关系为v’=2*exp(20-t/10)*exp(-20*exp(-t/10));
t1=0:0、1:100;
v1=2*exp(20-20-t1/10)、*exp(-20*exp(-t1/10)); plot(t1,v1)
6、(生态系统的振荡现象)第一次世界大战中,因为战争很少捕鱼,按理战后应能捕到更多的鱼才就是。
可就是大战后,在地中海却捕不到鲨鱼,因而渔民大惑不解。
令x1为鱼饵的数量,x2为鲨鱼的数量,t为时间。
常微分方程组为
式中a1, a2, b1, b2都就是正常数。
第一式鱼饵x1的增长速度大体上与x1成正比,即按a1x1比率增加, 而被鲨鱼吃掉的部分按b1x1x2的比率减少;第二式中鲨鱼的增长速度由于生存竞争的自然死亡或互相咬食按a2x2的比率减少,但又根据鱼饵的量的变化按b1x1x2的比率增加。
对a1=3, b1=2, a2=2、5, b2=1, x1(0)=x2(0)=1求解。
画出解曲线图与相轨线图,可以观
察
到鱼饵与鲨鱼数量的周期振荡现象。
代入a1=3, b1=2, a2=2、5, b2=1, x1(0)=x2(0)=1,x1’=3x1-2x1x2, x2’=-2、5x2+x1x2;
function Dx=fun(t,x)
Dx=zeros(2,1);
Dx(1)=x(1)-2*x(1)*x(2);
Dx(2)=-2、5*x(2)+x(1)*x(2);
f=f(:);
[t,x]=ode15s('fun',[0,10],[1,1]);
plot(t,x);。