PAGE\*MERGEFORMAT14开课学院、实验室：DS1421实验时间:2013年06月01日课程名称计算方法实验项目名称常微分方程的数值解法实验项目类型验证演示综合设计其他指导教师胡小兵成绩√1．实验目的：（1）学会四阶龙格-库塔方法的使用（2）设计出相应的算法，编制相应的函数子程序（3）会用这些函数解决实际问题2．实验内容（1）分别取h=0.05，N=10；h=0.025，N=20；h=0.01，N=50，用四阶龙格-库塔方法求解微分方程初值问题：y’=-50y，y(0)=10（2）某跳伞者在t=0时刻从飞机上跳出，假设初始时刻的垂直速度为0，且跳伞者垂直下落。已知空气阻力为F=cv2，其中c为常数，v为垂直速度，向下方方向为正。写出此跳伞者的速度满足的微分方程；若此跳伞者的质量为M=70kg，且已知c=0.27kg/m，利用四阶龙格-库塔公式计算t<=20s的速度（取h=0.1s）3.实验结果及分析四阶龙格-库塔方法%四阶经典R-K公式作数值计算clc;clearall;F='-50*y';a=0;b=0.5;%先对h=0.025计算h=0.025;%步长n=(b-a)/h;X1=a:h:b;Y1=zeros(1,n+1);Y1(1)=10;fori=1:nx=X1(i);y=Y1(i);K1=h*eval(F);x=x+h/2;y=y+K1/2;K2=h*eval(F);x=x;y=Y1(i)+K2/2;K3=h*eval(F);x=X1(i)+h;y=Y1(i)+K3;K4=h*eval(F);Y1(i+1)=Y1(i)+(K1+2*K2+2*K3+K4)/6;end%再对h=0.01计算h=0.01;%步长n=(b-a)/h;X2=a:h:b;Y2=zeros(1,n+1);Y2(1)=10;fori=1:nx=X2(i);y=Y2(i);K1=h*eval(F);x=x+h/2;y=y+K1/2;K2=h*eval(F);x=x;y=Y2(i)+K2/2;K3=h*eval(F);x=X2(i)+h;y=Y2(i)+K3;K4=h*eval(F);Y2(i+1)=Y2(i)+(K1+2*K2+2*K3+K4)/6;end%再对h=0.05计算h=0.05;%步长n=(b-a)/h;X3=a:h:b;Y3=zeros(1,n+1);Y3(1)=10;fori=1:nx=X3(i);y=Y3(i);K1=h*eval(F);x=x+h/2;y=y+K1/2;K2=h*eval(F);x=x;y=Y3(i)+K2/2;K3=h*eval(F);x=X3(i)+h;y=Y3(i)+K3;K4=h*eval(F);Y3(i+1)=Y3(i)+(K1+2*K2+2*K3+K4)/6;end%准确解temp=[];f=dsolve('Dy=-50*y','y(0)=10','x');df=zeros(1,n+1);fori=1:n+1temp=subs(f,'x',X2(i));df(i)=double(vpa(temp));endX1disp('h=0.025');disp(Y1);X2disp('h=0.01');disp(Y2);X3disp('h=0.05');disp(Y3);disp('准确值');df程序运行得到如下结果:X1=Columns1through1000.02500.05000.07500.10000.12500.15000.17500.20000.2250Columns11through200.25000.27500.30000.32500.35000.37500.40000.42500.45000.4750Column210.5000h=0.025Columns1through1010.00003.07450.94530.29060.08940.02750.00840.00260.00080.0002Columns11through200.00010.00000.00000.00000.00000.00000.00000.00000.00000.0000Column210.0000X2=Columns1through1000.01000.02000.03000.04000.05000.06000.07000.08000.0900Columns11through200.10000.11000.12