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

Matlab 应用 使用Euler和Rungkutta方法解臂状摆的能量方程 背景 单摆是常见的物理模型，为了得到摆角θ的关于时间的函数，来描述单摆运动。由角动量定理我们知道 化简得到 在一般的应用和计算中，只考虑摆角在5度以内的小摆动，因为可以吧简化为，这样比较容易解。实际上这是一个解二阶常微分方程的问题。 在这里的单摆是一种特别的单摆，具有均匀的质量M分布在长为2的臂状摆上， 使用能量法建立方程 化简得到 重力加速度取9.80665 1使用欧拉法 令，这样降阶就把二阶常微分方程转化为一阶微分方程组，再利用向前Euler方法数值求解。 y(i+1)=y(i)+h*z(i); z(i+1)=z(i)+h*7.35499*cos(y(i)); y(0)=0 z(0)=0 精度随着h的减小而更高，因为向前欧拉方法的整体截断误差与h同阶，（因为是用了泰勒公式）所以欧拉方法的稳定区域并不大。 RK4-四阶龙格库塔方法 使用四级四阶经典显式Rungkutta公式 稳定性很好，RK4法是四阶方法，每步的误差是h5阶，而总积累误差为h4阶。所以比欧拉稳定。 运行第三个程序：在一幅图中显示欧拉法和RK4法，随着截断误差的积累，欧拉法产生了较大的误差 h=0.01 h=0.0001，仍然是开始较为稳定，逐渐误差变大 总结：RK4是很好的方法，很稳定，而且四阶是很常用的方法，因为到五阶的时候精度并没有相应提升。通过这两种方法计算出角度峰值y=3.141593，周期是1.777510。 三个程序 欧拉法 clear; clc h=0.00001; a=0;b=25; x=a:h:b; y(1)=0; z(1)=0; for i=1:length(x)-1 % 欧拉 y(i+1)=y(i)+h*z(i); z(i+1)=z(i)+h*7.35499*cos(y(i)); end plot(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(欧拉); 龙格库塔方法 先定义函数rightf_sys1.m function w=rightf_sys1(x,y,z) w=7.35499*cos(y); clear; clc; %set(0,RecursionLimit,500) h=0.01; a=0;b=25; x=a:h:b; RK_y(1)=0; %初值 RK_z(1)=0; %初值 for i=1:length(x)-1 K1=RK_z(i); L1=rightf_sys1(x(i),RK_y(i),RK_z(i)); % K1 and L1 K2=RK_z(i)+0.5*h*L1; L2=rightf_sys1(x(i)+0.5*h,RK_y(i)+0.5*h*K1,RK_z(i)+0.5*h*L1); K3=RK_z(i)+0.5*h*L2; L3=rightf_sys1(x(i)+0.5*h,RK_y(i)+0.5*h*K2,RK_z(i)+0.5*h*L2); K4=RK_z(i)+h*L3; L4=rightf_sys1(x(i)+h,RK_y(i)+h*K3,RK_z(i)+h*L3); % K4 and L4 RK_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); end plot(x,RK_y,b+); xlabel(Variable x); ylabel(Variable y); A=[x,RK_y]; [y,T]=max(RK_y); legend(RK4方法); fprintf(角度峰值等于%d,y) %角度的峰值也就是π fprintf(\n)