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)
- 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);。