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

Matlab 应用

使用Euler和Rungkutta方法解臂状摆的能量方程

背景 单摆是常见的物理模型，为了得到摆角θ的关于时间的函数，来描述单

摆运动。由角动量定理我们知道

M=Jε

化简得到

dθg

+sinθ=02

dtl

在一般的应用和计算中，只考虑摆角在5度以内的小摆动，因为可以吧简化为θ，

这样比较容易解。实际上这是一个解二阶常微分方程的问题。

在这里的单摆是一种特别的单摆，具有均匀的质量M分布在长为2的臂状摆上， 使用能量法建立方程

2

T=W

化简得到

12

Jω=mg∆h

2

dθ

=7.35499cosθ2dt

重力加速度取9.80665

2

1使用欧拉法

dy令z=，这样降阶就把二阶常微分方程转化为一阶微分方程组，再利用向前

dx

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));

%龙格库塔

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); % K2 and L2

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); % K3 and L3

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,Euler_y,'r-',x,RK_y,'b-'); [y,T]=max(RK_y);

fprintf('角度峰值等于%d',y) %角度的峰值也就是π fprintf('\n')

fprintf('周期等于%d',T*h) %周期 xlabel('时间'); ylabel('角度'); legend('欧拉','RK4');

Matlab 应用

使用Euler和Rungkutta方法解臂状摆的能量方程

背景 单摆是常见的物理模型，为了得到摆角θ的关于时间的函数，来描述单

摆运动。由角动量定理我们知道

M=Jε

化简得到

dθg

+sinθ=02

dtl

在一般的应用和计算中，只考虑摆角在5度以内的小摆动，因为可以吧简化为θ，

这样比较容易解。实际上这是一个解二阶常微分方程的问题。

在这里的单摆是一种特别的单摆，具有均匀的质量M分布在长为2的臂状摆上， 使用能量法建立方程

2

T=W

化简得到

12

Jω=mg∆h

2

dθ

=7.35499cosθ2dt

重力加速度取9.80665

2

1使用欧拉法

dy令z=，这样降阶就把二阶常微分方程转化为一阶微分方程组，再利用向前

dx

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));

%龙格库塔

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); % K2 and L2

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); % K3 and L3

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,Euler_y,'r-',x,RK_y,'b-'); [y,T]=max(RK_y);

fprintf('角度峰值等于%d',y) %角度的峰值也就是π fprintf('\n')

fprintf('周期等于%d',T*h) %周期 xlabel('时间'); ylabel('角度'); legend('欧拉','RK4');


相关文章

  • 数值分析习题大作业

    • 一. [设计题一] .............................................................................................................. ...查看


  • 欧拉法与龙格库塔法比较分析

    • 解微分方程的欧拉法,龙格-库塔法简单实例比较 欧拉方法(Euler method)用以对给定初值的常微分方程(即初值问题)求解分为前EULER法.后退EULER法.改进的EULER法. 缺点: 欧拉法简单地取切线的端点作为下一步的起点进行计 ...查看


  • 微分方程求解数值方法

    • 重庆三峡学院数值计算方法实验课测验 微分方程求解的后退欧拉法(后退梯形法).龙格库塔 法(三阶.四阶公式) 院 系 数学与统计学院 专 业 信息与计算科学 姓 名 庞欢 年 级 2009级 学 号 [1**********]2 指导教师 刘 ...查看


  • 算法大全第15章常微分方程的解法

    • 第十五章 常微分方程的解法 建立微分方程只是解决问题的第一步,通常需要求出方程的解来说明实际现象,并加以检验.如果能得到解析形式的解固然是便于分析和应用的,但是我们知道,只有线性常系数微分方程,并且自由项是某些特殊类型的函数时,才可以肯定得 ...查看


  • 多种微分方程数值计算方法分析

    • 2010年8月 第4期 城 市 勘 测 Urban Geotechnical Investigation& Surveying Aug. 2010 No. 4 文章编号:1672-8262(2010)04-117-03中图分类号:O ...查看


  • 实验四 常微分方程初值问题数值解法

    • 实 验 报 告 课程名称 数值分析 实验项目 常微分方程问题初值问题数值解法 一. 实验目的 1.理解如何在计算机上实现用Euler 法.改进Euler 法.Runge -Kutta 算法求一阶常微分方程初值问题 ⎧y '(x ) =f ( ...查看


  • 悬崖跳水水池深度问题

    • 悬崖跳水水池深度问题 摘 要 本文探讨悬崖跳水水池深度设定问题,以实现水池深度设定既保证运动员人身安全且使成本消耗最低为目标.对此问题,将建立物理模型,运用物理学.流体力学知识,结合微分方程进行求解. 对于问题一,本文将建立物理模型,将跳水 ...查看


  • 电力系统稳定计算软件的开发

    • 电力系统稳定计算软件的开发 摘要 随着世界各国的电力系统发展的越来越庞大,电力系统的运行也随之越来越复杂,系统中发生的事故越来越难以用传统的分析方法预测,对电力系统仿真技术提出了更高的要求.目前,国际上有多种电力系统分析软件,在各国的电力系 ...查看


  • 6常微分方程的求解

    • §6常微分方程的求解 一.知识背景 常微分方程在物理学,工程技术中运用非常广泛,相当重要.许多物质运动的过程用常微分方程来描述,如:质点的加速运动.谐振动.电容充放电过程及电感通电断电等过程,因此,求解常微分方程成为很多物理问题求解的一种常 ...查看


热门内容