数学实验基础 实验报告1常微分方程Word下载.docx
《数学实验基础 实验报告1常微分方程Word下载.docx》由会员分享,可在线阅读,更多相关《数学实验基础 实验报告1常微分方程Word下载.docx(11页珍藏版)》请在冰豆网上搜索。
[x,y]
ans=
01.00000.10001.10000.20001.22000.30001.3620
0.40001.52820.50001.72100.60001.94310.70002.1974
0.80002.48720.90002.81591.00003.18751.10003.6062
1.20004.07691.30004.60451.40005.19501.50005.8545
1.60006.58991.70007.40891.80008.31981.90009.3318
2.000010.45502.100011.70052.200013.08052.300014.6086
2.400016.29952.500018.16942.600020.23642.700022.5200
2.800025.04202.900027.82623.000030.8988
ode45:
[x,y]=ode45('
[0,3],1)
01.00000.05021.05280.10051.11090.15071.1746
0.20101.24420.27601.35960.35101.48990.42601.6361
0.50101.79960.57601.98170.65102.18380.72602.4074
0.80102.65440.87602.92640.95103.22541.02603.5535
1.10103.91311.17604.30651.25104.73641.32605.2056
1.40105.71721.47606.27441.55106.88101.62607.5406
1.70108.25741.77609.03591.85109.88081.926010.7974
2.001011.79122.076012.86832.151014.03512.226015.2986
2.301016.66642.376018.14662.451019.74782.526021.4796
2.601023.35222.676025.37642.751027.56412.826029.9281
2.901032.48202.925733.36942.950534.27962.975235.2134
3.000036.1711
解析解:
y=dsolve('
Dy=x+y'
'
y(0)=1'
x'
)
y=
2*exp(x)-x-1
(2)
functionf=Fun(t,y)
f=[y
(2);
0.01*y
(2)^2-2*y
(1)+sin(t)];
[t,y]=euler('
[0,5],[0,1],0.2)
[t,y]=ode45('
[0,5],[0,1])
t=
00.00010.00010.00020.00020.00050.00070.00100.00120.0025
0.00370.00500.00620.01250.01880.02510.03130.06270.09410.1255
0.15690.28190.40690.53190.65690.78190.90691.03191.15691.2819
1.40691.53191.65691.78191.90692.03192.15692.28192.40692.5319
2.65692.78192.90693.03193.15693.28193.40693.53193.65693.7819
3.90694.03194.15694.28194.40694.53194.65694.74274.82854.9142
5.0000
y=
01.00000.00011.00000.00011.00000.00021.00000.00021.0000
0.00051.00000.00071.00000.00101.00000.00121.00000.00251.0000
0.00371.00000.00501.00000.00621.00000.01251.00000.01881.0000
0.02510.99990.03130.99980.06270.99870.09410.99650.12530.9934
0.15640.98930.27860.96320.39660.92200.50850.86620.61260.7967
0.70720.71460.79080.62100.86200.51760.91980.40580.96320.2876
0.99150.16471.00430.03921.0013-0.08690.9826-0.21170.9485-0.3331
0.8996-0.44900.8365-0.55780.7605-0.65770.6725-0.74710.5742-0.8246
0.4669-0.88890.3525-0.93930.2327-0.97480.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.99450.0300-0.98830.1146-0.97480.1985
-0.95430.2809
2.求一通过原点的曲线,它在
处的切线斜率等于
若
上限增为1.58,1.60会发生什么?
f=2*x+y.^2;
[0,1.57],0)
x=
00.03930.07850.11780.15700.19630.23550.27480.31400.35330.39250.43180.47100.51030.54950.58880.62800.6673
0.70650.74580.78500.82430.86350.90280.94200.98131.02051.05981.09901.13831.17751.21681.25601.29531.33451.37381.41301.42481.43671.44851.46041.47221.48401.49591.50771.51401.52031.52651.53281.53761.54241.54721.55191.55431.55671.55911.56141.56311.56471.56641.56811.56851.56901.56951.5700
00.00150.00620.01390.02470.03860.05560.07580.09920.12590.15590.18950.22660.26750.31240.36150.41520.47380.53780.60760.68410.76790.86010.96201.07511.20141.34341.50451.68921.90372.15572.45772.82823.30033.90564.73175.95496.44317.01167.68328.49029.482110.717012.309014.455115.922017.708019.939022.816425.645029.228233.967340.591044.943450.308857.122966.108774.310884.712398.4901117.7875124.9206132.9699142.1268152.6415
上限增为1.58,1.60,则超出运算的范围,发生溢出。
3.求解刚性方程组:
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('
[0,50],[1,-1])
t=00.00000.00000.00010.00010.00020.00030.00040.00050.00060.00060.00070.00080.00090.00100.00120.00130.00140.00150.00160.00180.00190.00200.00210.00230.00240.00250.00270.00280.00300.00310.00320.00330.00350.00360.00370.00390.00400.00420.00430.00450.00460.00490.00520.00550.00570.00600.00630.00680.00730.00790.00840.00890.00940.01070.01190.01320.01440.01570.01690.02160.02620.03080.03550.04010.06870.09730.12600.15460.33260.51060.68860.86651.04451.51461.98472.45482.92493.39493.86504.70505.54516.38517.22518.06518.905110.253311.601612.949814.298015.646216.994520.309423.624326.939330.254233.569136.884141.884146.884150.0000
1.0000-1.00000.9653-0.96530.9317-0.93170.8994-0.89930.7486-0.74840.6225-0.62220.5175-0.51720.4462-0.44570.3847-0.38420.3316-0.33110.2860-0.28530.2314-0.23070.1874-0.18650.1518-0.15080.1229-0.12190.0996-0.09850.0784-0.07710.0617-0.06030.0487-0.04720.0384-0.03680.0304-0.02870.0239-0.02210.0189-0.01690.0150-0.01290.0120-0.00970.0096-0.00720.0078-0.00530.0062-0.00350.0050-0.00220.0041-0.00120.0036-0.00050.0032-0.00000.00290.00040.00270.00070.00260.00100.00250.00120.00240.00150.00230.00170.00230.00190.00230.00200.00240.00210.00240.00220.00250.00240.00260.00260.00270.00270.00290.00280.00300.00300.00310.00310.00340.00340.00370.00370.00390.00390.00420.00420.00440.00440.00470.00470.00530.00530.00600.00600.00660.00660.00720.00720.00780.00780.00840.00840.01070.01070.01300.01300.01530.01530.01760.01760.01990.01990.03380.03380.04750.04750.06100.06100.07440.07440.15320.15320.22530.22530.29130.29130.35160.35160.40680.40680.53110.53110.62930.62930.70700.70700.76840.76840.81690.81690.85530.85530.90510.90510.93780.93780.95930.95930.97330.97330.98250.98250.98850.98850.99430.99430.99730.99730.99870.99870.99930.99930.99960.99960.99980.99981.00031.00031.00041.00041.00011.00010.99990.99991.00001.00001.00001.00001.00011.00011.00001.00001.00001.0000
4.(温度过程)夏天把开有空调的室内一支读数为20℃的温度计放到户外,10分钟后读25.2℃,再过10分钟后读数28.32℃。
建立一个较合理的模型来推算户外温度。
解:
设
并且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'
fun=
Inlinefunction:
fun(c,t)=c
(1)+exp(-c
(2)*t)*(-c
(1)+20)
lsqcurvefit(fun,[301],[1020],[25.228.32])
Localminimumfound.
Optimizationcompletedbecausethesizeofthegradientislessthan
thedefaultvalueofthefunctiontolerance.
<
stoppingcriteriadetails>
33.00000.0511
得户外温度c=33,比例系数k=0.0511
5.(广告效应)某公司生产一种耐用消费品,市场占有率为5%时开始做广告,一段时间的市场跟踪调查后,该公司发现:
单位时间内购买人口百分比的相对增长率与当时还没有买的百分比成正比,且估得此比例系数为0.5。
(1)建立该问题的数学模型,并求其数值解与模拟结果作以比较;
设单位时间购买人口百分比为:
x,则有数学模型:
则解析解:
模拟结果为:
functionf=Fun(t,x)
f=0.5*(1-x).*x;
[t,x]=ode45('
[0,10],0.05)
t=
00.10580.21150.31730.42310.67310.92311.17311.42311.67311.92312.17312.4231
2.67312.92313.17313.42313.67313.92314.17314.42314.67314.92315.17315.4231
5.67315.92316.17316.42316.67316.92317.17317.42317.67317.92318.17318.4231
8.67318.92319.17319.42319.56739.71159.855810.0000
0.05000.05260.05530.05810.06110.06860.07710.08640.09680.10830.12100.1349
0.15020.16690.18500.20460.22570.24830.27230.29780.32460.35250.38160.4115
0.44200.47310.50430.53550.56640.59680.62650.65520.68290.70930.73440.7581
0.78020.80090.82010.83780.85410.86290.87120.87900.8865
plot(x,t,'
.'
holdon
fplot(inline('
2*(log(x)-log(1-x)-log(0.05)+log(0.95))'
),[0.05,0.9])
holdoff
(2)厂家问:
要做多少时间广告,可使市场购买率达到80%?
n=min(find(x>
0.8));
t(n)
8.6731
6.(肿瘤生长)肿瘤大小V生长的速率与V的a次方成正比,其中a为形状参数,0a1;
而其比例系数K随时间减小,减小速率又与当时的K值成正比,比例系数为环境参数b。
设某肿瘤参数a=1,b=0.1,K的初始值为2,V的初始值为1。
问
(1)此肿瘤生长不会超过多大?
建立模型:
functionf=Fun(t,v)
a=1;
b=0.1;
f
(1)=v
(2).*v
(1).^a;
f
(2)=b.*v
(2);
[t,v]=ode45('
[0,50],[1,2]);
[t,k]=ode45('
plot(t,v,'
r-'
plot(t,k,'
b.'
max(v)
Inf296.8273
(2)过多长时间肿瘤大小翻一倍?
(3)何时肿瘤生长速率由递增转为递减?
(4)若参数a=2/3呢?
a=2/3;
选做题:
1.(生态系统的振荡现象)第一次世界大战中,因为战争很少捕鱼,按理战后应能捕到更多的鱼才是。
可是大战后,在地中海却捕不到鲨鱼,因而渔民大惑不解。
令x1为鱼饵的数量,x2为鲨鱼的数量,t为时间。
微分方程为
(5.20)
式中a1,a2,b1,b2都是正常数。
第一式鱼饵x1的增长速度大体上与x1成正比,即按a1x1比率增加,而被鲨鱼吃掉的部分按b1x1x2的比率减少;
第二式中鲨鱼的增长速度由于生存竞争的自然死亡或互相咬食按a2x2的比率减少,但又根据鱼饵的量的变化按b2x1x2的比率增加。
对a1=3,b1=2,a2=2.5,b2=1,x1(0)=x2(0)=1求解。
画出解曲线图和相轨线图,可以观察到鱼饵和鲨鱼数量的周期振荡现象。
2.编写四阶Runge-Kutta法程序。
functionvarargout=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('
ix(i)y(i)\n'
n
fprintf('
%2d%4.4f%4.4f\n'
i,x(i),y(i));
functionz=f(x,y)
z=-2*x*y^2;
function[y,x]=lgkt4j(x0,xn,y0,h)
x=x0:
xn;
n=length(x);
y1=x;
y1
(1)=y0;
fori=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;