数学实验基础 实验报告(1)常微分方程

实验一 常微分方程
1. 分别用Euler 法和ode45解下列常微分方程并与解析解比较： (1) ,(0)1,13y x y y x '=+=<<
Euler 法：
function [t,y]=euler(Fun,tspan,y0,h) t=tspan(1):h:tspan(2); y(1)=y0;
for i=1:length(t)-1
y(i+1)=y(i)+h.*feval(Fun,t(i),y(i)); end t=t'; y=y';
function f=Fun(x,y) % 常微分方程的右端函数 f=x+y;
>> [x,y]=euler('Fun',[0,3],1,0.1)
>> [x,y] ans =
0 1.0000 0.1000 1.1000 0.2000 1.2200 0.3000 1.3620 0.4000 1.5282 0.5000 1.7210 0.6000 1.9431 0.7000 2.1974 0.8000 2.4872 0.9000 2.8159 1.0000 3.1875 1.1000 3.6062 1.2000 4.0769 1.3000 4.6045 1.4000 5.1950 1.5000 5.8545 1.6000 6.5899 1.7000 7.4089 1.8000 8.3198 1.9000 9.3318 2.0000 10.4550 2.1000 11.7005 2.2000 13.0805 2.3000 14.6086 2.4000 16.2995 2.5000 18.1694 2.6000 20.2364 2.7000 22.5200 2.8000 25.0420 2.9000 27.8262 3.0000 30.8988
ode45：
>> [x,y]=ode45('Fun',[0,3],1) ans =
0 1.0000 0.0502 1.0528 0.1005 1.1109 0.1507 1.1746
0.2010 1.2442 0.2760 1.3596 0.3510 1.4899 0.4260 1.6361
0.5010 1.7996 0.5760 1.9817 0.6510 2.1838 0.7260 2.4074
实验一 常微分方程
0.8010 2.6544 0.8760 2.9264 0.9510 3.2254 1.0260 3.5535
1.1010 3.9131 1.1760 4.3065 1.2510 4.7364 1.3260 5.2056
1.4010 5.7172 1.4760 6.2744 1.5510 6.8810 1.6260 7.5406
1.7010 8.2574 1.7760 9.0359 1.8510 9.8808 1.9260 10.7974
2.0010 11.7912 2.0760 12.8683 2.1510 14.0351 2.2260 15.2986
2.3010 16.6664 2.3760 18.1466 2.4510 19.7478 2.5260 21.4796
2.6010 2
3.3522 2.6760 25.3764 2.7510 27.5641 2.8260 29.9281
2.9010 32.4820 2.9257 3
3.3694 2.9505 3
4.2796 2.9752 3
5.2134
3.0000 36.1711
解析解：>> y=dsolve('Dy=x+y','y(0)=1','x') y =
2*exp(x) - x - 1
(2) 2
0.01()2sin(),(0)0,(0)1,05y y y t y y t ''''-+===<< Euler 法：
实验一常微分方程
function f=Fun(t,y)
% 常微分方程的右端函数
f=[y(2);0.01*y(2)^2-2*y(1)+sin(t)];
>> [t,y]=euler('Fun',[0,5],[0,1],0.2)
ode45：
>> [t,y]=ode45('Fun',[0,5],[0,1])
t =
0 0.0001 0.0001 0.0002 0.0002 0.0005 0.0007 0.0010 0.0012 0.0025
0.0037 0.0050 0.0062 0.0125 0.0188 0.0251 0.0313 0.0627 0.0941 0.1255
0.1569 0.2819 0.4069 0.5319 0.6569 0.7819 0.9069 1.0319 1.1569 1.2819
1.4069 1.5319 1.6569 1.7819 1.9069
2.0319 2.1569 2.2819 2.4069 2.5319
2.6569 2.7819 2.9069
3.0319 3.1569 3.2819 3.4069 3.5319 3.6569 3.7819
3.9069
4.0319 4.1569 4.2819 4.4069 4.5319 4.6569 4.7427 4.8285 4.9142
5.0000
y =
0 1.0000 0.0001 1.0000 0.0001 1.0000 0.0002 1.0000 0.0002 1.0000
0.0005 1.0000 0.0007 1.0000 0.0010 1.0000 0.0012 1.0000 0.0025 1.0000
0.0037 1.0000 0.0050 1.0000 0.0062 1.0000 0.0125 1.0000 0.0188 1.0000
0.0251 0.9999 0.0313 0.9998 0.0627 0.9987 0.0941 0.9965 0.1253 0.9934
0.1564 0.9893 0.2786 0.9632 0.3966 0.9220 0.5085 0.8662 0.6126 0.7967
0.7072 0.7146 0.7908 0.6210 0.8620 0.5176 0.9198 0.4058 0.9632 0.2876
0.9915 0.1647 1.0043 0.0392 1.0013 -0.0869 0.9826 -0.2117 0.9485 -0.3331
0.8996 -0.4490 0.8365 -0.5578 0.7605 -0.6577 0.6725 -0.7471 0.5742 -0.8246
实验一 常微分方程
0.4669 -0.8889 0.3525 -0.9393 0.2327 -0.9748 0.1095 -0.9950 -0.0154 -0.9996
-0.1398 -0.9887 -0.2619 -0.9624 -0.3798 -0.9212 -0.4916 -0.8657 -0.5957 -0.7970
-0.6904 -0.7161 -0.7742 -0.6242 -0.8460 -0.5228 -0.9046 -0.4134 -0.9491 -0.2978
-0.9789 -0.1777 -0.9934 -0.0549 -0.9945 0.0300 -0.9883 0.1146 -0.9748 0.1985
-0.9543 0.2809
2. 求一通过原点的曲线，它在(,)x y 处的切线斜率等于2
2,0 1.57.x y x +<<若x 上限增为1.58，1.60会发生什么？
function f=Fun(x,y) % 常微分方程的右端函数 f=2*x+y.^2;
>> [x,y]=ode45('Fun',[0,1.57],0) x =
0 0.0393 0.0785 0.1178 0.1570 0.1963 0.2355 0.2748 0.3140 0.3533 0.3925 0.4318 0.4710 0.5103 0.5495 0.5888 0.6280 0.6673 0.7065 0.7458 0.7850 0.8243 0.8635 0.9028 0.9420 0.9813 1.0205 1.0598 1.0990 1.1383 1.1775 1.2168 1.2560 1.2953 1.3345 1.3738 1.4130 1.4248 1.4367 1.4485 1.4604 1.4722 1.4840 1.4959 1.5077 1.5140 1.5203 1.5265 1.5328 1.5376 1.5424 1.5472 1.5519 1.5543 1.5567 1.5591 1.5614 1.5631 1.5647 1.5664 1.5681 1.5685 1.5690 1.5695 1.5700 y =
实验一 常微分方程
0 0.0015 0.0062 0.0139 0.0247 0.0386 0.0556 0.0758 0.0992
0.1259 0.1559 0.1895 0.2266 0.2675 0.3124 0.3615 0.4152 0.4738 0.5378 0.6076 0.6841 0.7679 0.8601 0.9620 1.0751 1.2014 1.3434 1.5045 1.6892 1.9037 2.1557 2.4577 2.8282 3.3003 3.9056 4.7317 5.9549 6.4431 7.0116 7.6832 8.4902 9.4821 10.7170 12.3090 14.4551 15.9220 17.7080 19.9390 22.8164 25.6450 29.2282 33.9673 40.5910 44.9434 50.3088 57.1229 66.1087 74.3108 84.7123 98.4901 117.7875 124.9206 132.9699 142.1268 152.6415
00.20.40.60.81 1.2 1.4 1.6
若x 上限增为1.58，1.60,则超出运算的范围，发生溢出。

3. 求解刚性方程组：
11212
121100.25999.750.5,(0)1,050.
999.751000.250.5,(0)1,y y y y x y y y y '=-++=⎧<<⎨
'=-+=-⎩
function f=Fun(t,y)
% 常微分方程的右端函数
f(1)=-1000.25*y(1)+999.75*y(2)+0.5; f(2)=999.75*y(1)-1000.25*y(2)+0.5; f=f(:);
[t,y]=ode15s('Fun',[0,50],[1,-1])
t = 0 0.0000 0.0000 0.0001 0.0001 0.0002 0.0003 0.0004 0.0005 0.0006 0.0006
0.0007 0.0008 0.0009 0.0010 0.0012 0.0013 0.0014 0.0015 0.0016 0.0018 0.0019 0.0020 0.0021 0.0023 0.0024 0.0025 0.0027 0.0028 0.0030 0.0031 0.0032 0.0033 0.0035 0.0036 0.0037 0.0039 0.0040 0.0042 0.0043 0.0045 0.0046 0.0049 0.0052
实验一常微分方程
0.0055 0.0057 0.0060 0.0063 0.0068 0.0073 0.0079 0.0084 0.0089 0.0094 0.0107
0.0119 0.0132 0.0144 0.0157 0.0169 0.0216 0.0262 0.0308 0.0355 0.0401 0.0687
0.0973 0.1260 0.1546 0.3326 0.5106 0.6886 0.8665 1.0445 1.5146 1.9847 2.4548
2.9249
3.3949 3.8650
4.7050
5.5451
6.3851
7.2251
8.0651 8.9051 10.2533 11.6016
12.9498 14.2980 15.6462 16.9945 20.3094 23.6243 26.9393 30.2542 33.5691 36.8841 41.8841
46.8841 50.0000
y =
1.0000 -1.0000 0.9653 -0.9653 0.9317 -0.9317 0.8994 -0.8993 0.7486 -0.7484 0.6225
-0.6222 0.5175 -0.5172 0.4462 -0.4457 0.3847 -0.3842 0.3316 -0.3311 0.2860 -0.2853
0.2314 -0.2307 0.1874 -0.1865 0.1518 -0.1508 0.1229 -0.1219 0.0996 -0.0985 0.0784
-0.0771 0.0617 -0.0603 0.0487 -0.0472 0.0384 -0.0368 0.0304 -0.0287 0.0239 -0.0221
0.0189 -0.0169 0.0150 -0.0129 0.0120 -0.0097 0.0096 -0.0072 0.0078 -0.0053 0.0062
-0.0035 0.0050 -0.0022 0.0041 -0.0012 0.0036 -0.0005 0.0032 -0.0000 0.0029 0.0004
0.0027 0.0007 0.0026 0.0010 0.0025 0.0012 0.0024 0.0015 0.0023 0.0017 0.0023
0.0019 0.0023 0.0020 0.0024 0.0021 0.0024 0.0022 0.0025 0.0024 0.0026 0.0026
0.0027 0.0027 0.0029 0.0028 0.0030 0.0030 0.0031 0.0031 0.0034 0.0034 0.0037
0.0037 0.0039 0.0039 0.0042 0.0042 0.0044 0.0044 0.0047 0.0047 0.0053 0.0053
0.0060 0.0060 0.0066 0.0066 0.0072 0.0072 0.0078 0.0078 0.0084 0.0084 0.0107
0.0107 0.0130 0.0130 0.0153 0.0153 0.0176 0.0176 0.0199 0.0199 0.0338 0.0338
0.0475 0.0475 0.0610 0.0610 0.0744 0.0744 0.1532 0.1532 0.2253 0.2253 0.2913
0.2913 0.3516 0.3516 0.4068 0.4068 0.5311 0.5311 0.6293 0.6293 0.7070 0.7070
0.7684 0.7684 0.8169 0.8169 0.8553 0.8553 0.9051 0.9051 0.9378 0.9378 0.9593
0.9593 0.9733 0.9733 0.9825 0.9825 0.9885 0.9885 0.9943 0.9943 0.9973 0.9973
0.9987 0.9987 0.9993 0.9993 0.9996 0.9996 0.9998 0.9998 1.0003 1.0003 1.0004
1.0004 1.0001 1.0001 0.9999 0.9999 1.0000 1.0000 1.0000 1.0000 1.0001 1.0001
1.0000 1.0000 1.0000 1.0000
实验一 常微分方程
4. (温度过程)夏天把开有空调的室内一支读数为20℃的温度计放到户外，10分钟后读2
5.2℃, 再过10分钟后读数28.32℃。

建立一个较合理的模型来推算户外温度。

解：设)('T c k T -=并且T(0)=20, T(10)=25.2, T(20)=28.52,其中0<t<60.
>> dsolve('DT=k*(c-T)','T(0)=20','t') ans =
c - (c - 20)/exp(k*t)
利用T(10)=25.2, T(20)=28.52拟合： fun=inline('c(1)+exp(-c(2)*t)*(-c(1)+20)','c','t') fun =
Inline function:
fun(c,t) = c(1)+exp(-c(2)*t)*(-c(1)+20) >> lsqcurvefit(fun,[30 1],[10 20],[25.2 28.32]) Local minimum found.
Optimization completed because the size of the gradient is less than the default value of the function tolerance. <stopping criteria details> ans =
33.0000 0.0511
得户外温度c=33,比例系数k=0.0511
5. (广告效应)某公司生产一种耐用消费品，市场占有率为5%时开始做广告，一段时间的市场跟踪调查后，该公司发现：单位时间内购买人口百分比的相对增长率与当时还没有买的百分比成正比，且估得此比例系数为0.5。

实验一 常微分方程
(1) 建立该问题的数学模型，并求其数值解与模拟结果作以比较；
设单位时间购买人口百分比为：x ，则有数学模型：⎪⎩⎪⎨⎧=-=%
5)0()
1(*5.0*x dt dx
x x
则解析解：
)95.0ln 05.0ln )1ln((ln *2t +---=x x
模拟结果为：
function f=Fun(t,x) % 常微分方程的右端函数
f=0.5*(1-x).*x;
>> [t,x]=ode45('Fun',[0,10],0.05)
t =
0 0.1058 0.2115 0.3173 0.4231 0.6731 0.9231 1.1731 1.4231 1.6731 1.9231 2.1731 2.4231
2.6731 2.9231
3.1731 3.4231 3.6731 3.9231
4.1731 4.4231 4.6731 4.9231
5.1731 5.4231
5.6731 5.9231
6.1731 6.4231 6.6731 6.9231
7.1731 7.4231 7.6731 7.9231
8.1731 8.4231
8.6731 8.9231 9.1731 9.4231 9.5673 9.7115 9.8558 10.0000
x =
0.0500 0.0526 0.0553 0.0581 0.0611 0.0686 0.0771 0.0864 0.0968 0.1083 0.1210 0.1349
0.1502 0.1669 0.1850 0.2046 0.2257 0.2483 0.2723 0.2978 0.3246 0.3525 0.3816 0.4115
0.4420 0.4731 0.5043 0.5355 0.5664 0.5968 0.6265 0.6552 0.6829 0.7093 0.7344 0.7581
0.7802 0.8009 0.8201 0.8378 0.8541 0.8629 0.8712 0.8790 0.8865
>> plot(x,t,'.'); >> hold on
>> fplot(inline('2*(log(x)-log(1-x)-log(0.05)+log(0.95))'),[0.05,0.9])
实验一 常微分方程
>> hold off
(2) 厂家问：要做多少时间广告，可使市场购买率达到80%？ >> n=min(find(x>0.8)); >> t(n) ans = 8.6731
6. (肿瘤生长) 肿瘤大小V 生长的速率与V 的a 次方成正比，其中a 为形状参数，0≤a ≤1;而其比例系数K 随时间减小，减小速率又与当时的K 值成正比，比例系数为环境参数b 。

设某肿瘤参数a=1, b=0.1, K 的初始值为2，V 的初始值为1。

问 (1)此肿瘤生长不会超过多大？
建立模型：⎪⎪⎩⎪⎪⎨⎧====2)0(,1)0(,k bk dt
dk v kv dt
dv a
function f=Fun(t,v) % 常微分方程的右端函数
a=1; b=0.1;
f(1)=v(2).*v(1).^a;
实验一常微分方程
f(2)=b.*v(2);
f=f(:);
>> [t,v]=ode45('Fun',[0,50],[1,2]);
>> [t,k]=ode45('Fun',[0,50],[1,2]);
>> plot(t,v,'r-');
>> hold on
>> plot(t,k,'b.')
>> max(v)
ans =
Inf 296.8273
(2)过多长时间肿瘤大小翻一倍？
(3)何时肿瘤生长速率由递增转为递减？
(4)若参数a=2/3呢？
function f=Fun(t,v)
% 常微分方程的右端函数
a=2/3;
b=0.1;
f(1)=v(2).*v(1).^a;
f(2)=b.*v(2);
f=f(:);
>> [t,v]=ode45('Fun',[0,50],[1,2]);
>> [t,k]=ode45('Fun',[0,50],[1,2]);
>> plot(t,v,'r-');
实验一 常微分方程
11 / 12
2013春 数学实验 实验一 常微分方程
>> hold on
>> plot(t,k,'b.')
8
选做题：
1.(生态系统的振荡现象)第一次世界大战中，因为战争很少捕鱼，按理战后应能捕到更多的鱼才是。

可是大战后，在地中海却捕不到鲨鱼，因而渔民大惑不解。

令x 1为鱼饵的数量，x 2为鲨鱼的数量，t 为时间。

微分方程为
(5.20) 式中a 1, a 2, b 1, b 2都是正常数。

第一式鱼饵x 1的增长速度大体上与x 1成正比，即按a 1x 1比率增加, 而被鲨鱼吃掉的部分按b 1x 1x 2的比率减少；第二式中鲨鱼的增长速度由于生存竞争的自然死亡或互相咬食按a 2x 2的比率减少，但又根据鱼饵的量的变化按b 2x 1x 2的比率增加。

对a 1=3, b 1=2, a 2=2.5, b 2=1, x 1(0)=x 2(0)=1求解。

画出解曲线图和相轨线图，可以观察到鱼饵和鲨鱼数量的周期振荡现象。

2.编写四阶Runge-Kutta 法程序。

function varargout=saxplaxliu(varargin)
实验一常微分方程
clc;
clear
x0=0;
xn=1.2;
y0=1;
h=0.1;
[y,x]=lgkt4j(x0,xn,y0,h);
n=length(x);
fprintf(' i x(i) y(i)\n');
for i=1:n
fprintf('%2d %4.4f %4.4f\n',i,x(i),y(i));
end
function z=f(x,y)
z=-2*x*y^2;
function [y,x]=lgkt4j(x0,xn,y0,h)
x=x0:h:xn;
n=length(x);
y1=x;
y1(1)=y0;
for i=1:n-1
K1=f(x(i),y1(i));
K2=f(x(i)+h/2,y1(i)+h/2*K1);
K3=f(x(i)+h/2,y1(i)+h/2*K2);
K4=f(x(i)+h,y1(i)+h*K3);
y1(i+1)=y1(i)+h/6*(K1+2*K2+2*K3+K4);
end
y=y1;
12 / 12 2013春数学实验实验一常微分方程。