4阶Runge-Kutta法求解一阶常微分方程(可打印修改)

参考教材《数值分析》李乃成.梅立泉clearclcformat longm=input('请输入常微分方程的阶数m=');a=input('请输入x下限a=');b=input('请输入x上限b=');h=input('请输入步长h=');ym=input('令y(1,1)=y,y(2,1)=y’,y(3,1)=y’’...请输入ym=','s'); %输入的时候必须按照这个形式输入y1=y(1,1);if m==1 %一阶初值问题单独求解mm=(b-a)/h;y(1,1)=input('请输入在初值点的函数值f(a)=');x=a;y11(1)=y(1,1);for k1=2:(mm+1)y1=y(1,1);K(1,1)=h*(eval(ym)); %计算K1x=x+h/2;y(1,1)=y1+K(1,1)/2;y1=y(1,1);K(1,2)=h*(eval(ym)); %计算K2x=x;y(1,1)=y1+K(1,2)/2-K(1,1)/2;y1=y(1,1);K(1,3)=h*(eval(ym)); %计算K3x=x+h/2;y(1,1)=y1+K(1,3)-K(1,2)/2;y1=y(1,1);K(1,4)=h*(eval(ym)); %计算K4y11(k1)=y11(k1-1)+(K(1,1)+2*K(1,2)+2*K(1,3)+K(1,4))/6; y(1,1)=y11(k1);x=a+(k1-1)*h;endy11else %高阶初值问题mm=(b-a)/h; %一共要求解mm个数据点for k2=1:m %读取初值条件fprintf('请输入%d阶导数的初值f(%d)(a)=\n',(k2-1),(k2-1));y(k2,1)=input('=');endfor k2=1:my22(1,k2)=y(k2,1); %先把初值保存在矩阵y22(m,n)中,m表示第几个所求点,n表示第n阶初值endx=a;for k4=2:(mm+1) %求解mm个数据点的循环for k=1:(m-1) %计算K1,包括每一阶的K1 K(k,1)=h*y(k+1,1); %y(k+1,1)中k+1表示第k+1阶,1表示第一个点;K(k,1)中k表示阶数,1表示K1endK(m,1)=h*(eval(ym));x=x+h/2; %求解K1之前,先重新对x和y赋值for k3=1:my(k3,1)=y(k3,1)+K(k3,1)/2;endfor k=1:(m-1) %计算K2K(k,2)=h*y(k+1,1);endK(m,2)=h*(eval(ym));x=x;for k3=1:my(k3,1)=y(k3,1)-K(k3,1)/2+K(k3,2)/2;endfor k=1:(m-1) %计算K3K(k,3)=h*y(k+1,1);endK(m,3)=h*(eval(ym));x=x+h/2;for k3=1:my(k3,1)=y(k3,1)+K(k3,3)-K(k3,2)/2; %这里容易出错endfor k=1:(m-1) %计算K4K(k,4)=h*y(k+1,1);endK(m,4)=h*(eval(ym));for k5=1:my22(k4,k5)=y22(k4-1,k5)+(K(k5,1)+2*K(k5,2)+2*K(k5,3)+K(k5,4))/6; %这里,除了要求出下一个点的数值,还要求出相应的导数值endfor k6=1:m %除了对y(1,1)重新赋值外,还要对y(2,1)等重新赋值y(k6,1)=y22(k4,k6);endx=a+(k4-1)*h;endy22(:,1) end。

计算方法上机作业——龙格库塔法龙格库塔法（Runge-Kutta method）是一种常用于求解常微分方程（Ordinary Differential Equation，ODE）的数值解法。

它是由德国数学家卡尔·龙格（Carl Runge）和马丁·威尔海姆·库塔（Martin Wilhelm Kutta）在20世纪初提出的。

龙格库塔法的基本思想是通过数值逼近来计算微分方程的近似解。

在讲解龙格库塔法之前，我们先来简单回顾一下ODE的一阶常微分方程的基本形式：y′(y)=y(y,y)其中，y(y,y)是已知函数。

龙格库塔法的核心是使用差分逼近计算函数的斜率。

假设我们要求解的方程为：y′(y)=y(y,y),y(y)=y₀所需计算的点为y₀,y₁,...,yy，对应的函数值为y₀,y₁,...,yy，其中y是步长的个数。

龙格库塔法通过递推关系式来计算估计值，并不断更新当前点的函数值。

接下来以龙格库塔法的经典四阶形式为例进行说明。

该方法的基本方程如下：yy+1=yy+(y₁+2y₂+2y₃+y₄)/6y₁=ℎy(yy,yy)y₂=ℎy(yy+ℎ/2,yy+y₁/2)y₃=ℎy(yy+ℎ/2,yy+y₂/2)y₄=ℎy(yy+ℎ,yy+y₃)其中y表示当前步骤，ℎ表示步长，yy表示当前点的函数值，y₁,y₂,y₃和y₄则表示对应的斜率。

使用龙格库塔法，我们可以通过不断递归计算来求得指定区间（例如[y,y]）上的函数值。

具体步骤如下：1.确定求解区间[y,y]和初始点(y₀,y₀)以及步长ℎ。

2.初始化：设置yy=y₀，yy=y₀。

3.对所有y=0,...,y−1：计算y₁,y₂,y₃和y₄，根据上述递推关系式。

根据递推关系式计算yy+1更新当前点的函数值，即yy+1=y(yy+1)。

更新当前点的y值，即yy+1=yy+ℎ。

4.返回结果：最终求得的函数值。

需要注意的是，选择适当的步长对最终结果的精度和计算效率都有重要影响。

《数值分析》课程设计常微分方程初值问题数值解的实现和分析—四阶Runge-kutta方法及预估-校正算法常微分方程初值问题数值解的实现和分析—四阶Runge-kutta方法及预估-校正算法摘要求解常微分方程的初值问题，Euler方法，改进的Euler方法及梯形方法精度比较低，所以本文构造高精度单步的四级Runge-kutta方法及高精度的多步预估—校正算法及其Matlab编程来实现对常微分方程初值问题的求解，使在求解常微分方程时，对以前积分方法的收敛速度及精度都有了很高的提高。

关键词：Runge-kutta方法，Adams方法，预估—校正算法，Matlab目录1.前言 (1)2. 几个简单的数值积分法 (2)2.1Runge-kutta方法 (2)2.1.1 Runge-kutta方法的应用 (5)2.2预估—校正算法 (7)2.2.1 Adams数值积分方法简介及预估—校正算法 (7)2.2.2 预估—校正算法的应用 (12)3. 结果分析 (16)总结 (17)参考文献 (18)英文原文和中文翻译 (19)1英文原文 (19)2中文翻译 (20)1.前言常微分方程的初值问题是微分方程定解问题的一个典型代表，以下面的例子介绍常微分方程初值问题数值解的基本思想和原理。

例1.1 一重量垂直作用于弹簧所引起的震荡，当运动阻力与速度的平方成正比时，可借助如下二阶常微分方程描述若令和，则上述二阶常微分方程可化成等价的一阶常微分方程组类似于例1.1，对于m阶常微分方程其中。

若定义可得如下等价的一阶常微分方程组我们知道多数常微分方程主要靠数值解法。

所谓数值解法，就是寻求解在一系列离散节点上的近似值。

相邻两个节点之间的间距称为步长[1]。

2. 几个简单的数值积分法2.1 Runge-kutta方法Runge在1985年提出了一种基于Euler折线法的新的数值方法，此后这种新的数值方法又经过其同胞K.Heun和Kutta的努力[2]，发展完善成为后世所称的Runge-kutta 方法。

matlab经典的4级4阶runge kutta法-回复使用MATLAB 实现经典的4 阶4 级Runge-Kutta 法引言：数值计算是现代科学和工程中的一个重要领域，它涉及到通过计算机模拟来解决数学问题。

在数值计算中，求解微分方程是一个常见的任务。

Runge-Kutta 法是求解微分方程的一种常见方法，它可以用于数值求解常微分方程和偏微分方程。

本文将介绍经典的4 级4 阶Runge-Kutta 法的原理，并使用MATLAB 来实现该方法。

一、原理介绍：Runge-Kutta 法是数值计算领域中最常用的方法之一。

它通过将微分方程的解逐步逼近来求解微分方程。

经典的4 级4 阶Runge-Kutta 法基于以下公式：\begin{align*}k_1 &= h f(t_n, y_n) \\k_2 &= h f(t_n + \frac{h}{2}, y_n + \frac{k_1}{2}) \\k_3 &= h f(t_n + \frac{h}{2}, y_n + \frac{k_2}{2}) \\k_4 &= h f(t_n + h, y_n + k_3) \\y_{n+1} &= y_n + \frac{1}{6}(k_1 + 2k_2 + 2k_3 + k_4)\end{align*}其中，h 是步长，t_n 是当前时间点，y_n 是当前的解，f(t, y) 是微分方程的右手函数。

二、算法实现：现在我们将使用MATLAB 实现经典的4 级4 阶Runge-Kutta 法，并解决一个简单的一阶常微分方程。

首先，我们定义一个MATLAB 函数，用于实现4 级4 阶Runge-Kutta 法。

函数接受输入参数为微分方程的右手函数f(t, y)，初始时间t_0，初始解y_0，以及步长h。

函数输出为一个数组，包含了每个时间点的解。

以下是MATLAB 代码实现：matlabfunction y = runge_kutta(f, t0, y0, h, num_steps)初始化解数组y = zeros(num_steps+1, 1);y(1) = y0;循环计算每个时间点的解for i = 1:num_stepst = t0 + (i-1)*h;计算k1, k2, k3, 和k4k1 = h * f(t, y(i));k2 = h * f(t + h/2, y(i) + k1/2);k3 = h * f(t + h/2, y(i) + k2/2);k4 = h * f(t + h, y(i) + k3);计算下一个时间点的解y(i+1) = y(i) + (k1 + 2*k2 + 2*k3 + k4)/6;endend接下来，我们使用这个函数来解决一个简单的一阶常微分方程。