数值分析Matlab作业龙格库塔欧拉方法解二阶微分方程
- 格式:doc
- 大小:201.50 KB
- 文档页数:12
Matlab 应用
使用Euler 和Rungkutta 方法解臂状摆的能量方程
背景 单摆是常见的物理模型,为了得到摆角θ的关于时间的函数,来描述单摆运动。由角动量定理我们知道
ε
J M =
化简得到 0sin 22
=+θθl
g dt d 在一般的应用和计算中,只考虑摆角在5度以内的小摆动,因为可以吧简化为θ,这样比较容易解。实际上这是一个解二阶常微分方程的问题。
在这里的单摆是一种特别的单摆,具有均匀的质量M 分布在长为2的臂状摆上, 使用能量法建立方程
W
T = h mg ∆=2J 2
1ω 化简得到 θθcos 35499.722=dt
d 重力加速度取9.80665
1使用欧拉法 令dx
dy z =,这样降阶就把二阶常微分方程转化为一阶微分方程组,再利用向前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 同阶,(因为是用了泰勒公式)所以欧拉方法的稳定区域并不大。
2.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')
fprintf('周期等于%d',T*h) %周期
两个方法在一起对比
使用跟上一个相同的函数rightf_sys1.m
clear;
clc; %清屏
h=0.0001;
a=0;b=25;
x=a:h:b;
Euler_y(1)=0;
Euler_z(1)=0; %欧拉的初值RK_y(1)=0;
RK_z(1)=0; %龙格库塔初值for i=1:length(x)-1
%先是欧拉法
Euler_y(i+1)=Euler_y(i)+h*Euler_z(i);
Euler_z(i+1)=Euler_z(i)+h*7.35499*cos(Euler_y(i));
%龙格库塔