当前位置：文档之家 › 四阶龙格库塔实验报告

四阶龙格库塔实验报告

  • 格式：docx
  • 大小：64.14 KB
  • 文档页数：4

三、四阶Runge-Kutta 法求解常微分方程
一、龙格库塔法的思想
根据第九章的知识可知道，Euler 方法的局部截断误差是2()O h ，而当用Euler 方法估计出
1,()
(1)n n n n y y hf x y +=+ 再用梯形公式
111[(,)(,)](2)2n n n n n n h y y f x y f x y +++=++
进行校正，即采用改进Euler 方法得出数值解的截断误差为3()O h 。

由Lagrange 微分中值定理
'11()()()()()(,())(3)n n n n n y x y x y x x y x hf y ξξξ++=+-=+ 记*(,())k hf y ξξ=，得到
*1()()(4)n n y x y x k +=+
这样只要给出一种计算*k 的算法，就能得到相应的计算公式。

用这种观点的来分析Euler 方法和改进Euler 方法，Euler 方法的迭代公式可改写为
11
1(,)
n n n n y y k k hf x y +=+=
改进Euler 方法的预报-校正公式可改写为 1121211()2
(,),(,)
n n n n n n y y k k k hf x y k hf x h y k +=++==++ Euler 方法实际上是用一个点处的值1k 近似*k ，而改进Euler 方法是用两个点处的值1k ，和2k ，做算术平均值近似*k 自然改进Euler 方法要优于Euler 方法。

因此，可以想到假如在1[,]n n x x +内多预报几个点值i k ，并用他们的加权平均值作为*k 的近似值，则有可能构造出具有更高精度的计算公式，这就是Runge-Kutta 法的基本思想。

二、四阶龙格库塔法
由Runge-Kutta 的基本思想，构造四阶Runge-Kutta 法是利用1234,,k k k k 和的加权平均值来近似*k ，因此令
111223344121113
2211243312213(,)(,)(5)(,)(,)n n n n n n n n n n y y w K w K w K w K K hf x y K hf x h y K K hf x h y K K K hf x h y K K K αβαβγαβγη+=++++⎧⎪=⎪⎪=++⎨⎪=+++⎪⎪=++++⎩
使得
511()()n n y x y O h ++-=
即其总体截断误差为4()O h 。

采用泰勒公式展开，经过复杂的推导，得到一个具有13个参数，11个方程的线性方程组。

由于方程的个数少于未知量的个数，因此方程有无穷多个解。

可以根据情况得到几种常用的解，即得到相应的四阶公式。

最常见的四阶公式如式（6）： ，112341213243(22)6(,)
1(,)(6)221(,)22(,)
n n n n n n n n n n h y y K K K K K hf x y h K hf x y K h K hf x y K K hf x h y K +⎧=++++⎪⎪=⎪⎪⎪=++⎨⎪⎪=++⎪⎪=++⎪⎩
也称为标准四阶Runge-Kutta 法。

三、四阶龙格库塔法程序说明及应用
龙格库塔的计算程序
function [x,y] =Runge(ydot_fun,x0,y0,h,N )
x=zeros(1,N+1);y=zeros(length(y0),N+1);
x(1)=x0;y(:,1)=y0;
for n=1:N
x(n+1)=x(n)+h;
k1=h*feval(ydot_fun,x(n),y(:,n));
k2=h*feval(ydot_fun,x(n)+1/2*h,y(:,n)+1/2*k1);
k3=h*feval(ydot_fun,x(n)+1/2*h,y(:,n)+1/2*k2);
k4=h*feval(ydot_fun,x(n)+h,y(:,n)+k3);
y(:,n+1)=y(:,n)+1/6*(k1+2*k2+2*k3+k4);
end
程序解释及使用
该算法可以对一阶微分方程，一阶微分方程组进行有效的求解。

ydot_fun 为一阶微分方程的函数，x0为初始点，y0为初始向量，h 为步长，N 为区间的等分数，x 为Xn 构成的向量，y 为Yn 构成的矩阵。

程序调用方法：
1，先编写要求解的一阶微分方程或方程组的函数文件文件，将该文件和Runge 文件放到同一个目录下。

2. 调用求解程序，[x,y]=Runge(@dot_fun,x0,y0,h,N)，运行后即可得出结果。

或者用内部函数调用：
输入：ydot_fun=(x,y)[]
[x,y]= Runge(ydot_fun,x0,y0,h,N)
实例求解
课本304页题目：
用标准4级4阶R-K 法求解
，'''''''''23(0)1,(0)3,(0)2y y y y x y y y ⎧=+-+-⎪
⎨=-==⎪⎩
， 取步长h=，计算(1)y 的近似值，并与解析解()21x y x xe x =+-作比较。

解：首先将三阶方程改写成微分方程组的形式：
令'''''
12132,,y y y y y y y y =====得如下微分方程组 '12'23'''3211'''1
2323(0)1,(0)3,(0)2y y y y y y y y x y y y ⎧=⎪=⎪⎨=+-+-⎪⎪=-==⎩ 在文件中编写待求解微分方程组，调用计算程序，保留5位小数得：
表3-1 三阶微分方程求解结果
表3-1中第二行1y 为原三阶微分方程对应的数值解，第二行2y 为其一阶导数值，第三行3y 为其二阶导数值。

由结果可知，数值解y(1)=，其对应的精确解析解(1)1y e =+=3.的相对误差为，可知四阶龙格库塔法具有很高的代数精度。

合集下载

Runge-Kutta报告

Runge-Kutta报告

实验报告——四级四阶Runge-Kutta 公式1110500158 刘雪婷一 实验问题从上一个报告中知Euler 法和梯形法是最简单的求解常微分方程离散点值的方法。

但这两种方法的精确度都比较低。

在解微分方程的时候为了使精度更高，我们使用了Taylor 级数法，通过增加展开项数来提高方法的阶。

但是，这种方法计算量大，所以我们用一种间接使用Taylor 展开式来获取高阶的方法。

这种方法就是Runge-Kutta 法，既可以减少计算量，又可以提高精度。

二 实验方法{ du dt =ℱ(t,u)u (t 0)=u 0* 对*式左右两边积分，得到： u (t n+1)−u (t n )=∫ℱ(t,u(t))dt t n+1t n 积分中值定理1 ⇒ h ℱ(t n +θh,u(t n +θh))我们用［t n ,t n+1］中若干个点的线性组合来近似θ。

u n+1=u (t n )+h ∑ωi k i s i=1#k i =ℱ(t n +hαi ,u (t n )+h ∑βi,j k j i−1j=1) ,i ≥1 k 1=ℱ(t n ,u(t n )) 将#式两边在t n 处Taylor 展开，然后按照h 的幂排列，让同幂的两项系数相等，以此来确定ωi αj βij 。

利用上述方法，我们得到了古典Rung-Kutta 公式（四级四阶）： y n+1=y n +h 6(K 1+2K 2+2K 3+K 4)K 1=ℱ(x n ,y n )K 2=ℱ(x n +h 2,y n +h 2K 1 K 3=ℱ(x n +h 2,y n +h 2K 2 K 4=ℱ(x n +h,y n +hK 3 三 初始量{dy dx =x 2−y y (0)=10≤x ≤1 h=0.1四 结论可以看到，RK 方法得到的结果精确到小数点后5位。

五 结论分析当h=0.2时，结果精确到小数点后4位。

当h=0.01时，结果精确到了小数点后9位。

龙格—库塔 C++

龙格—库塔 C++

实验四 龙格-库塔法【实验内容】1. 用标准四阶龙格－库塔方法设计算法;2. 编程解微分方程初值问题.3.在课后习题中选择一个题目编程计算。

交回实验报告与计算结果。

【实验方法或步骤】四阶龙格－库塔方法的计算步骤求解'0()(,)()()y x f x y a x b y a y ⎧=≤≤⎨=⎩ 对上述给定的(,)f x y ，用四阶龙格-库塔法求解常微分方程初值问题112341213243(22)6(,)11(,)2211(,)22(,)n n n n n n n n n n h y y k k k k k hf x y k hf x h y k k hf x h y k k hf x h y k +⎧=++++⎪⎪=⎪⎪⎪=++⎨⎪⎪=++⎪⎪=++⎪⎩[实验程序]（自编程序 .cpp 文件）#include<stdio.h>double f(double a, double b){double c;c = a + b;return c;}int main(void)inti;double a, b, m, n, h;double k1, k2, k3, k4;printf("请输入端点值a和b：\n");scanf("%lf %lf", &a, &b);printf("请输入整数n：\n");scanf("%lf", &n);printf("请输入初始值：\n");scanf("%lf", &m);h = (b-a)/n;printf("在区间（%3.1lf %3.1lf）上步长为%3.1lf初值为: \n", a, b, h);printf("X0 = %lf, Y0 = %lf\n", a ,m);for (i=0; i<n; i++){k1 = f(a, m);k2 = f(a+h/2, m+h*k1/2);k3 = f(a+h/2, m+h*k2/2);k4 = f(a+h, m+h*k3);m = m + (k1+2*k2+2*k3+k4)*h/6;a = a + h;printf("X%d = %lf, Y%d = %lf\n",i+1, a ,i+1, m);}return 0;【实验结果】（课后4题）。

龙格库塔实验报告

龙格库塔实验报告

一、实验背景常微分方程（ODE）在自然科学、工程技术等领域中具有广泛的应用。

然而，许多微分方程无法得到精确解析解，因此需要借助数值方法进行求解。

龙格-库塔（Runge-Kutta）方法是一种常用的数值求解常微分方程的方法，具有精度高、稳定性好等优点。

本实验旨在通过编写程序，实现四阶龙格-库塔方法，并验证其在求解常微分方程中的有效性和准确性。

二、实验目的1. 理解四阶龙格-库塔方法的基本原理和计算步骤。

2. 编写程序实现四阶龙格-库塔方法。

3. 选取典型常微分方程，验证四阶龙格-库塔方法的求解精度和稳定性。

三、实验原理四阶龙格-库塔方法是一种基于泰勒级数展开的数值方法，其基本思想是将微分方程的解在某个区间内进行近似，并通过迭代计算得到近似解。

具体步骤如下：1. 初始化：给定初始条件y0，步长h，求解区间[a, b]。

2. 迭代计算：对于k=1, 2, ..., n（n为迭代次数），- 计算k1 = f(xk-1, yk-1)（f为微分方程的右端函数）；- 计算k2 = f(xk-1 + h/2, yk-1 + h/2 k1)；- 计算k3 = f(xk-1 + h/2, yk-1 + h/2 k2)；- 计算k4 = f(xk-1 + h, yk-1 + h k3)；- 更新y值：yk = yk-1 + (h/6) (k1 + 2k2 + 2k3 + k4)；- 更新x值：xk = xk-1 + h；3. 输出结果：输出最终的近似解y(n)。

四、实验步骤1. 编写程序实现四阶龙格-库塔方法。

2. 选取典型常微分方程，如：- y' = -y，初始条件y(0) = 1，求解区间[0, 2π]；- y' = y^2，初始条件y(0) = 1，求解区间[0, 1]。

3. 对每个常微分方程，设置不同的步长h和迭代次数n，分别计算近似解y(n)。

4. 将计算得到的近似解与解析解进行比较，分析四阶龙格-库塔方法的精度和稳定性。

四阶龙格-库塔(R-K)方法求常微分方程

四阶龙格-库塔(R-K)方法求常微分方程

中南大学MATLAB程序设计实践材料科学与工程学院2013年3月26日一、编程实现“四阶龙格－库塔（R-K）方法求常微分方程”，并举一例应用之。

【实例】采用龙格-库塔法求微分方程：,0)(1'00x x y y y 1、算法说明：在龙格-库塔法中，四阶龙格-库塔法的局部截断误差约为o(h5)，被广泛应用于解微分方程的初值问题。

其算法公式为：)22(63211k k k h y y nn其中：),()21,21()21,21(),(3423121hk y h x f k hk y h x f k hk y h x f k y x f k n nn n n n n n 2、流程图：2.1、四阶龙格－库塔（R-K ）方法流程图：2.2、实例求解流程图：输入待求微分方程、求解的自变量范围、初值以及求解范围内的取点数等。

确定求解范围内的步长k = 取点数？否求解：),()21,21()21,21(),(3423121hk y h x f k hk y h x f k hk y h x f k y x f k n nnn n n n n 求解并输出：)22(63211k k k h y y nn是结束算法开始输入求解的自变量范围求出待求简单微分方程的真值解用MA TLAB自带函数ode23求解待求微分方程用自编函数四阶龙格－库塔（R-K）方法求解待求微分方程结束3、源程序代码3.1、四阶龙格－库塔（R-K）方法源程序：function [x,y] = MyRunge_Kutta(fun,x0,xt,y0,PointNum,varargin)%Runge-Kutta 方法解微分方程形为 y'(t)=f(x,y(x))%此程序可解高阶的微分方程。

只要将其形式写为上述微分方程的向量形式%函数 f(x,y): fun%自变量的初值和终值：x0, xt%y0表示函数在x0处的值，输入初值为列向量形式%自变量在[x0,xt]上取的点数：PointNum%varargin为可输入项，可传适当参数给函数f(x,y)%x：所取的点的x值%y：对应点上的函数值if nargin<4 | PointNum<=0PointNum=100;endif nargin<3y0=0;endy(1,:)=y0(:)'; %初值存为行向量形式h=(xt-x0)/(PointNum-1); %计算步长x=x0+[0:(PointNum-1)]'*h; %得x向量值for k=1:(PointNum) %迭代计算f1=h*feval(fun,x(k),y(k,:),varargin{:});f1=f1(:)'; %得公式k1 f2=h*feval(fun,x(k)+h/2,y(k,:)+f1/2,varargin{:});f2=f2(:)'; %得公式k2 f3=h*feval(fun,x(k)+h/2,y(k,:)+f2/2,varargin{:});f3=f3(:)'; %得公式k3 f4=h*feval(fun,x(k)+h,y(k,:)+f3,varargin{:});f4=f4(:)'; %得公式k4 y(k+1,:)=y(k,:)+(f1+2*(f2+f3)+f4)/6; %得y(n+1) end3.2、实例求解源程序：%运行四阶R-K法clear, clc %清除内存中的变量x0=0;xt=2;Num=100;h=(xt-x0)/(Num-1);x=x0+[0:Num]*h;a=1;yt=1-exp(-a*x); %真值解fun=inline('-y+1','x','y'); %用inline构造函数f(x,y)y0=0; %设定函数初值PointNum=5; %设定取点数[x1,y1]=ode23(fun,[0,2],0);[xr,yr]=MyRunge_Kutta(fun,x0,xt,y0,PointNum);MyRunge_Kutta_x=xr'MyRunge_Kutta_y=yr'plot(x,yt,'k',x1,y1,'b--',xr,yr,'r-')legend('真值','ode23','Rung-Kutta法解',2)hold onplot(x1,y1,'bo',xr,yr,'r*')4、程序运行结果：MyRunge_Kutta_x =0 0.5000 1.0000 1.5000 2.0000MyRunge_Kutta_y =0 0.3932 0.6318 0.7766 0.8645二、变成解决以下科学计算问题：（一）[例7-2-4] 材料力学复杂应力状态的分析——Moore 圆。