计算方法实验报告-常微分方程的数值解法.docVIP

重 庆 大 学 学 生 实 验 报 告 实验课程名称 计算方法 开课实验室 DS1421 学 院 年级 专业 学 生 姓 名 学 号 开 课 时 间 至 学年 第 学期 总 成 绩 教师签名  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.实验过程 解： %RungeKutta.m %四阶经典R-K公式作数值计算 clc; F=-50*y; a=0; b=0.5; h=0.05; n=(b-a)/h; X=a:h:b; Y=zeros(1,n+1); Y(1)=10; for i=1:n x=X(i); y=Y(i); K1=h*eval(F); x=x+h/2; y=y+K1/2; K2=h*eval(F); x=x; y=Y(i)+K2/2; K3=h*eval(F); x=X(i)+h; y=Y(i)+K3; K4=h*eval(F); Y(i+1)=Y(i)+(K1+2*K2+2*K3+K4)/6; end %准确解 temp=[]; f=dsolve(Dy=-50*y,y(0)=10,x); df=zeros(1,n+1); for i=1:n+1 temp=subs(f,x,X(i)); df(i)=double(vpa(temp)); end disp(步长四阶经典R-K法准确值); disp([X,Y,df]); %画图观察效果 figure; plot(X,df,k*,X,Y,--r); grid on; title(四阶经典R-K法解常微分方程); legend(准确值,四阶经典R-K法); h=0.05，N=10，结果如下所示： (需将程序中，结果如下所示： (需将程序中，结果如下所示： 牛顿第二定律可得：70，c=0.27，初始条件v(0)=0 即： 所程序如下： %sudu.m clc; F=9.8-0.27/70*v^2; a=0; b=20; h=0.1; n=(b-a)/h; T=a:h:b; V=zeros(1,n+1); V(1)=0; for i=1:n t=T(i); v=V(i); K1=h*eval(F); t=t+h/2; v=v+K1/2; K2=h*eval(F); t=t; v=V(i)+K2/2; K3=h*eval(F); t=T(i)+h; v=V(i)+K3; K4=h*eval(F); V(i+1)=V(i)+(K1+2*K2+2*K3+K4)/6; end %准确解 temp=[]; f=dsolve(Dv=9.8-0.27/70*v^2,v(0)=0); df=zeros(1,n+1); for i=1:n+1 temp=subs(f,t,T(i)); df(i)=double(vpa(temp)); end disp(步长四阶经典R-K法准确值); disp([T,V,df]); %画图观察效果 figure; plot(T,df,k*,T,V,--r); grid on; title(速度v与时间t的函数曲线); legend(准确值,四阶经典R-K法); 结果如下所示： 0 0 0 0.1000 0.9799 0.9799 0.2000 1.9590 1.9590 0.3000 2.9367 2.9367 0.4000 3.9121 3.9121 0.5000 4.8846 4