数值分析经典例题1.y' = y , x [0,1] ,y (0) =1 , h = 0.1。
1求解析解。
2 Eular法3 R-K法○1解析法在MATLAB命令窗口执行clear>> x=0:0.1:1;>> y=exp(x);>> c=[y]'c =1.0000000000000001.1051709180756481.2214027581601701.3498588075760031.4918246976412701.6487212707001281.8221188003905092.0137527074704772.2255409284924682.4596031111569502.718281828459046○2Euler法在Matlab中建立M文件如下：function [x,y]=euler1(dyfun,xspan,y0,h)x=xspan(1):h:xspan(2);y(1)=y0;for n=1:length(x)-1y(n+1)=y(n)+h*feval(dyfun,x(n),y(n));endx=x';y=y'在MATLAB命令窗口执行clear>> dyfun=inline('y+0*x');>> [x,y]=euler1(dyfun,[0,1],1,0.1);>> [x,y]得到ans =0 1.0000000000000000.100000000000000 1.1000000000000000.200000000000000 1.2100000000000000.300000000000000 1.3310000000000000.400000000000000 1.4641000000000000.500000000000000 1.6105100000000000.600000000000000 1.7715610000000000.700000000000000 1.9487171000000000.800000000000000 2.1435888100000000.900000000000000 2.3579476910000001.0000000000000002.593742460100000○3R-K法（龙格-库塔法）在本题求解中，采用经典4阶龙格-库塔法首先在Matlab的M文件窗口对4阶龙格-库塔算法进行编程：function [x,y]=RungKutta41(dyfun,x0,y0,h,N)x=zeros(1,N+1);y=zeros(1,N+1);x(1)=x0;y(1)=y0;for n=1:Nx(n+1)=x(n)+h;k1=h*feval(dyfun,x(n),y(n));k2=h*feval(dyfun,x(n)+h/2,y(n)+1/2*k1);k3=h*feval(dyfun,x(n)+h/2,y(n)+1/2*k2);k4=h*feval(dyfun,x(n+1)+h,y(n)+k3);y(n+1)=y(n)+(k1+2*k2+2*k3+k4)/6;end在MATLAB命令窗口执行clear>> dyfun=inline('y','x','y');>> [x,y]=RungKutta41(dyfun,0,1,0.1,10);>> c=[x;y]'得到c =0 1.0000000000000000.100000000000000 1.1051708333333330.200000000000000 1.2214025708506950.300000000000000 1.3498584970625380.400000000000000 1.4918242400806860.500000000000000 1.6487206385968380.600000000000000 1.8221179620919330.700000000000000 2.0137516265967770.800000000000000 2.2255395632923150.900000000000000 2.4596014137800711.0000000000000002.718279744135166 ○4绘图'解析法','Euler法','R-K法' 绘制如下在MATLAB命令窗口执行clear>> x=0:0.1:1;>> y1=exp(x);>> dyfun=inline('y+0*x');>> [x,y2]=euler1(dyfun,[0,1],1,0.1);>> dyfun=inline('y','x','y');>> [x,y3]=RungKutta41(dyfun,0,1,0.1,10);>> plot(x,y1,'*')hold onplot(x,y2,'g','LineWidth',2)plot(x,y3,'b','LineWidth',2)legend('解析法','Euler法','R-K法')2.一个具有1400kg初始重量的小火箭，带有1040kg的燃料，点燃后垂直向上运动，火箭内的燃料以18kg/s的速率燃烧，提供31000N的推力。
试求燃料燃尽之前火箭的运动规律。
（高度h，速度v，加速度a随时间变化的情况）假定空气阻力D ∝v2, D=k*v2, k=0.3822 N*S2/m2。
解：本题利用Matlab求解过程如下：在MATLAB命令窗口执行>> clear>> a=0;>> v=0;>> h=0;>> aa(1:59)=0;>> vv(1:59)=0;>> hh(1:59)=0;>> for t=0:58aa(t+1)=a;vv(t+1)=v;hh(t+1)=h;l1=(9.8*18*t-0.3822*v^2+31000-1400*9.8)/( 1400-18*t);l2=(9.8*18*(t+0.5)-0.3822*(v+0.5*l1)^2+31000-1400*9.8)/(1400-18*(t+0.5)); l3=(9.8*18*(t+0.5)-0.3822*(v+0.5*l2)^2+31000-1400*9.8)/(1400-18*(t+0.5)); l4=(9.8*18*(t+1)-0.3822*(v+l3)^2+31000-1400*9.8)/( 1400-18*(t+1));h=h+v+(l1+l2+l3)/6;v=v+(l1+2*l2+2*l3+l4)/6;a=l1;end>> t=0:1:58;○1输入命令>> plot(t, aa,'r','LineWidth',2)得燃料燃尽前火箭飞行加速度变化曲线运行>> [t;aa]'得火箭飞行中加速度具体数据ans =0 01.000000000000000 12.3428571428571432.000000000000000 12.5882391353900853.000000000000000 12.7500563343637684.000000000000000 12.8224167027991375.000000000000000 12.8017614062130046.000000000000000 12.6871474793833887.000000000000000 12.4803824997476558.000000000000000 12.1859927919667659.000000000000000 11.81102349364022110.000000000000000 11.36468637522382411.000000000000000 10.85788687081145212.000000000000000 10.30267294687005013.000000000000000 9.71165369600690214.000000000000000 9.09743453012283315.000000000000000 8.47210933379155216.000000000000000 7.84683954413271617.000000000000000 7.23153787358350718.000000000000000 6.63466226065696419.000000000000000 6.06311518646046620.000000000000000 5.52223569746864621.000000000000000 5.01586667252432222.000000000000000 4.54647789380093923.000000000000000 4.11532581499936424.000000000000000 3.72263289340209025.000000000000000 3.36777228575659226.000000000000000 3.04944701534642927.000000000000000 2.76585595755594228.000000000000000 2.51484187625120129.000000000000000 2.29401912245700430.000000000000000 2.10088043498854631.000000000000000 1.93288358686951732.000000000000000 1.78751947040075033.000000000000000 1.66236369478680034.000000000000000 1.55511397232646135.000000000000000 1.46361557331463936.000000000000000 1.38587700395899337.000000000000000 1.32007785928313838.000000000000000 1.26457056388743739.000000000000000 1.21787746534880840.000000000000000 1.17868450607912341.000000000000000 1.14583248035056042.000000000000000 1.11830668929079343.000000000000000 1.09522563961762444.000000000000000 1.07582929101099345.000000000000000 1.05946724020834746.000000000000000 1.04558713435446347.000000000000000 1.03372352881725748.000000000000000 1.02348734264703549.000000000000000 1.01455601541566350.000000000000000 1.00666442996568051.000000000000000 0.99959663463942352.000000000000000 0.99317837421125853.000000000000000 0.98727041970092754.000000000000000 0.98176267247720355.000000000000000 0.97656900677018056.000000000000000 0.97162280629547257.000000000000000 0.96687314469179458.000000000000000 0.962281555544803 ○2输入命令plot(t, vv,'g','LineWidth',2)得燃料燃尽前火箭飞行速度变化曲线运行>> [t;vv]'得火箭飞行中速度具体数据ans =0 01 12.47221033291382 25.14859865803973 37.94247477466464 50.76239405614845 63.51464315701586 76.10593656299687 88.44616545539628 100.451030111789 112.04439422860810 123.16022350210111 133.7440082908212 143.75361596642313 153.15956644632514 161.94476842077215 170.10378919257316 177.64175495809617 184.57298992609918 190.91950284240219 196.70942048350120 201.97545237231521 206.75345230015122 211.0811********23 214.99689048355324 218.53896445078725 221.74457738991926 224.64939977519627 227.28711136337228 229.68911067400429 231.88434170507530 233.89921726903931 235.75761953723932 237.48096024095133 239.08828517202934 240.59640991767735 242.020*********36 243.37211856489437 244.66363889223138 245.90417600157539 247.10187370470240 248.26363998409441 249.39529677733742 250.50171889061843 251.58696134821244 252.65437492945345 253.70670998331346 254.74620886322247 255.77468750774548 256.79360682026449 257.80413458414750 258.80719869830651 259.80353253892152 260.79371325234953 261.7781937666954 262.75732927884255 263.73139893347756 264.70062336266657 265.66517870215258 266.625207644371 ○3输入命令>> plot(t, hh,'b','LineWidth',2)得燃料燃尽前火箭飞行高度时间变化曲线运行>> [t;hh]'得火箭飞行过程中高度具体数据ans =0 01 6.215658520407712 25.01258223286623 56.55209543058494 100.9062607015555 158.0543438047046 227.8818819853917 310.182********8 404.66236384119710 628.59189496914611 757.0903********12 895.88845811881713 1044.3963050495414 1202.0006618702415 1368.0771*******16 1542.0012728424917 1723.158482871218 1910.9524599680419 2104.8120991452720 2304.1968382571821 2508.6005113101622 2717.5538316988823 2930.6256637048224 3147.4232604062425 3367.5916517920926 3590.8123616153527 3816.8016185809628 4045.3082097125229 4276.1111035555730 4509.0169500850931 4743.8575440999832 4980.4873203651333 5218.7809322900334 5458.6309517211735 5699.9457154829236 5942.6473345075837 6186.6698735498838 6431.9577033482939 6678.4640224283240 6926.1495423027741 7174.9813273933742 7424.9317793819743 7675.9777547316244 7928.0998*******45 8181.2815187518346 8435.50898170547 8690.7702969745248 8947.0552*******49 9204.3547433037750 9462.6610118392951 9721.9669252538552 9982.266053650854 10505.820688321655 10769.065479724356 11033.281902896557 11298.465204156458 11564.6107884638。