数值分析Matlab作业龙格库塔欧拉方法解二阶微分方程VIP

数值分析Matlab作业龙格库塔欧拉方法解二阶微分方程数值分析Matlab作业龙格库塔欧拉方法解二阶微分方程数值分析Matlab作业龙格库塔欧拉方法解二阶微分方程数值分析Matlab作业龙格库塔欧拉方法解二阶微分方程编制仅供参考审核批准生效日期地址：电话：传真:邮编:Matlab应用使用Euler和Rungkutta方法解臂状摆的能量方程背景单摆是常见的物理模型，为了得到摆角θ的关于时间的函数，来描述单摆运动。由角动量定理我们知道化简得到在一般的应用和计算中，只考虑摆角在5度以内的小摆动，因为可以吧简化为，这样比较容易解。实际上这是一个解二阶常微分方程的问题。在这里的单摆是一种特别的单摆，具有均匀的质量M分布在长为2的臂状摆上，使用能量法建立方程化简得到重力加速度取1使用欧拉法令，这样降阶就把二阶常微分方程转化为一阶微分方程组，再利用向前Euler方法数值求解。y(i+1)=y(i)+h*z(i);z(i+1)=z(i)+h**cos(y(i));y(0)=0z(0)=0精度随着h的减小而更高，因为向前欧拉方法的整体截断误差与h同阶，（因为是用了泰勒公式）所以欧拉方法的稳定区域并不大。RK4-四阶龙格库塔方法使用四级四阶经典显式Rungkutta公式稳定性很好，RK4法是四阶方法，每步的误差是h5阶，而总积累误差为h4阶。所以比欧拉稳定。运行第三个程序：在一幅图中显示欧拉法和RK4法，随着截断误差的积累，欧拉法产生了较大的误差h=h=，仍然是开始较为稳定，逐渐误差变大总结：RK4是很好的方法，很稳定，而且四阶是很常用的方法，因为到五阶的时候精度并没有相应提升。通过这两种方法计算出角度峰值y=，周期是。三个程序欧拉法clear;clch=;a=0;b=25;x=a:h:b;y(1)=0;z(1)=0;fori=1:length(x)-1%欧拉y(i+1)=y(i)+h*z(i);z(i+1)=z(i)+h**cos(y(i));endplot(x,y,'r*');xlabel('时间');ylabel('角度');A=[x,y];%y(find(y==max(y)))%Num=(find(y==max(y)))[y,T]=max(y);fprintf('角度峰值等于%d',y)%角度的峰值也就是πfprintf('\n')fprintf('周期等于%d',T*h)%周期legend('欧拉');龙格库塔方法先定义函数functionw=rightf_sys1(x,y,z)w=*cos(y);clear;clc;%set(0,'RecursionLimit',500)h=;a=0;b=25;x=a:h:b;RK_y(1)=0;%初值RK_z(1)=0;%初值fori=1:length(x)-1K1=RK_z(i);L1=rightf_sys1(x(i),RK_y(i),RK_z(i));%K1andL1K2=RK_z(i)+*h*L1;L2=rightf_sys1(x(i)+*h,RK_y(i)+*h*K1,RK_z(i)+*h*L1);K3=RK_z(i)+*h*L2;L3=rightf_sys1(x(i)+*h,RK_y(i)+*h*K2,RK_z(i)+*h*L2);K4=RK_z(i)+h*L3;L4=rightf_sys1(x(i)+h,RK_y(i)+h*K3,RK_z(i)+h*L3);%K4andL4RK_y(i+1)=RK_y(i)+1/6*h*(K1+2*K2+2*K3+K4);RK_z(i+1)=RK_z(i)+1/6*h*(L1+2*L2+2*L3+L4);endplot(x,RK_y,'b+');xlabel('Variablex');ylabel('Variabley');A=[x,RK_y];[y,T]=max(RK_y);legend('RK4方法');fprintf('角度峰值等于%d',y)%角度的峰值也就是πfprintf('\n')fprintf('周期等于%d',T*h)%周期两个方法在一起对比使用跟上一个相同的函数clear;clc;%清屏h=;a=0;b=25;x=a:h:b;Euler_y(1)=0;Euler_z(1)=0;%欧拉的初值RK_y(1)=0;RK_z(1)=0;%龙格库塔初值fori=1:length(x)-1%先是欧拉法Euler_y(i+1)=Euler_y(i)+h*Euler_z(i);Euler_z(i+1)=Euler_z(i)+h**cos(Euler_y(i));%龙格库塔K1=RK_z(i);L1=rightf_sys1(x(i),RK_y(i),RK_z(i));%K1andL1K2=RK_z(i)+*h*L1;L2=rightf_sys1(x(i)+*h,RK_y(i)+*h*K1,RK_z(i)+*h*L1);%K2andL2K3=RK_z(i)+*h*L2;L3=rightf_sys1(x(i)+*h,RK_y(i)+*h*K2,RK_z(i)+*h*L2);%K3andL3K4=RK_z(i)+h*L3;L4=rightf_sys1(x(i)+h,RK_y(i)+h*K3,RK_z(i)+h*L3);%K4andL4RK_y(i+1)=RK_y(i)+1/6*h*(K1+2*K2+2*K3+K4);RK_z(i+1)=RK_z(i)+1/6*h*(L1+2*L2+2*L3+L4);endplot(x,Euler_y,'r-',x,RK_y,'b-');[y,T]=max(RK_y);