数值分析Matlab作业龙格库塔欧拉方法解二阶微分方程Word文档格式.doc
- 文档编号:960201
- 上传时间:2023-04-29
- 格式:DOC
- 页数:7
- 大小:207.75KB
数值分析Matlab作业龙格库塔欧拉方法解二阶微分方程Word文档格式.doc
《数值分析Matlab作业龙格库塔欧拉方法解二阶微分方程Word文档格式.doc》由会员分享,可在线阅读,更多相关《数值分析Matlab作业龙格库塔欧拉方法解二阶微分方程Word文档格式.doc(7页珍藏版)》请在冰点文库上搜索。
令,这样降阶就把二阶常微分方程转化为一阶微分方程组,再利用向前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;
fori=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)%角度的峰值也就是π
\n'
)
周期等于%d'
T*h)%周期
legend('
欧拉'
龙格库塔方法
先定义函数rightf_sys1.m
functionw=rightf_sys1(x,y,z)
w=7.35499*cos(y);
clc;
%set(0,'
RecursionLimit'
500)
h=0.01;
RK_y
(1)=0;
%初值
RK_z
(1)=0;
%初值
length(x)-1
K1=RK_z(i);
L1=rightf_sys1(x(i),RK_y(i),RK_z(i));
%K1andL1
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);
%K4andL4
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);
plot(x,RK_y,'
b+'
Variablex'
Variabley'
A=[x,RK_y];
[y,T]=max(RK_y);
RK4方法'
两个方法在一起对比
使用跟上一个相同的函数rightf_sys1.m
%清屏
h=0.0001;
Euler_y
(1)=0;
Euler_z
(1)=0;
%欧拉的初值
%龙格库塔初值
%先是欧拉法
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));
%龙格库塔
L1=rightf_sys1(x(i),RK_y(i),RK_z(i));
%K1andL1
L2=rightf_sys1(x(i)+0.5*h,RK_y(i)+0.5*h*K1,RK_z(i)+0.5*h*L1);
%K2andL2
%K3andL3
L4=rightf_sys1(x(i)+h,RK_y(i)+h*K3,RK_z(i)+h*L3);
plot(x,Euler_y,'
r-'
x,RK_y,'
b-'
y)%角度的峰值也就是π
'
RK4'
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值 分析 Matlab 作业 龙格库塔欧拉 方法 解二阶 微分方程