matlab编的4阶龙格库塔法解微分方程的程序

合集下载

四阶龙格-库塔法求解常微分方程的初值问题-matlab通用程序

参考教材《数值分析》李乃成.梅立泉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。

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

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

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

其算法公式为：)22(63211k k k hy y n n +++=+其中：⎪⎪⎪⎩⎪⎪⎪⎨⎧++=++=++==) ,()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 n n n n n n n 2、流程图：2.1、四阶龙格－库塔（R-K ）方法流程图：2.2、实例求解流程图：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 圆。

matlab经典的4级4阶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接下来，我们使用这个函数来解决一个简单的一阶常微分方程。

微分方程的数值解法matlab(四阶龙格—库塔法)

解析解： x x x1 3 2(((ttt))) 0 .0 8 1 1 2 P k 8 0siw n t) (2 .6 3 0 3 3 P k 0siw n t) (0 .2 12 2 2 P k 0siw n t)(
第一个质量的位移响应时程
Y (t)A(Y t)P(t)
(2)
Y (t)A(Y t)P(t)
3. Matlab 程序(主程序：ZCX)
t0;Y0;h;N;P0,w; %输入初始值、步长、迭代次数、初始激励力；
for i = 1 : N
t1 = t0 + h
P=[P0*sin(w*t0);0.0;0.0]
%输入t0时刻的外部激励力
Van der Pol方程
% 子程序（程序名： dYdt.m ） function Ydot = dYdt (t, Y) Ydot=[Y(2);-Y(2)*(Y(1)^2-1)-Y(1)];
或写为
function Ydot = dYdt (t, Y) Ydot=zeros(size(Y)); Ydot(1)=Y(2); Ydot(2)=-Y(2)*(Y(1).^2-1)-Y(1)];
Solver解算指令的使用格式
说明：
t0：初始时刻；tN：终点时刻 Y0：初值； tol：计算精度
[t, Y]=solver (‘ODE函数文件名’, t0, tN, Y0, tol);
ode45
输出宗量形式
y1 (t0 )
Y
y1
(t1
)
y
1
(t
2
)
y2 (t0 )
y
2
(
t1
)
y
2
(
t

matlab经典的4级4阶runge kutta法

MATLAB是一种用于算法开发、数据分析、可视化和数值计算的高级技术计算语言和交互式环境。

作为一个强大的工具，MATLAB提供了许多数值计算方法，其中4级4阶Runge-Kutta方法就是其中之一。

1. Runge-Kutta方法简介Runge-Kutta方法是求解常微分方程（ODE）的数值方法之一。

在MATLAB中，用户可以使用内置的ode45函数来调用4级4阶Runge-Kutta方法。

具体来说，4级4阶Runge-Kutta方法是一种单步迭代方法，通过在每个步骤中计算斜率来逐步逼近解析解。

它的优点是数值稳定性好，适用于多种类型的微分方程。

2. Runge-Kutta方法的公式4级4阶Runge-Kutta方法的一般形式如下：$$k_1 = hf(t_n, y_n)$$$$k_2 = hf(t_n + \frac{1}{2}h, y_n + \frac{1}{2}k_1)$$$$k_3 = hf(t_n + \frac{1}{2}h, y_n + \frac{1}{2}k_2)$$$$k_4 = hf(t_n + h, y_n + k_3)$$$$y_{n+1} = y_n + \frac{1}{6}(k_1 + 2k_2 + 2k_3 + k_4)$$其中，$t_n$是当前的独立变量值，$y_n$是当前的解向量，h是步长，$f(t_n, y_n)$是给定点$(t_n, y_n)$处的斜率。

通过不断迭代上述公式，可以逐步求解微分方程的数值解。

3. MATLAB中的4级4阶Runge-Kutta方法的应用在MATLAB中，用户可以使用ode45函数调用4级4阶Runge-Kutta方法来求解常微分方程。

使用ode45函数的基本语法如下：```matlab[t, y] = ode45(odefun, tspan, y0)```其中，odefun是用户定义的ODE函数句柄，tspan指定了求解的时间范围，y0是初始条件。

4阶经典龙格库塔公式求解微分方程

4阶经典龙格库塔公式求解微分方程4阶龙格库塔法（也称为四阶Runge-Kutta方法）是一种常用于求解常微分方程的数值方法。

它是由德国数学家卡尔·龙格以及他的学生马丁·威尔海姆·库塔于1901年独立提出的。

以下是详细的介绍。

1.问题描述我们考虑一个典型的常微分方程初值问题：dy/dx = f(x, y)，y(x0) = y0其中，x0是给定的初始点，y(x)是我们要求解的函数，f(x,y)是已知的函数。

2.方法原理四阶龙格库塔方法的基本思想是通过使用全区间的函数信息来估计下一步的函数值。

在每个步骤中，我们计算四个不同的估计量，然后将它们组合起来以获取最终的解。

具体来说，我们首先计算在初始点x0上的斜率k1=f(x0,y0)。

然后我们计算在x0+h/2处的斜率k2=f(x0+h/2,y0+h/2*k1)，其中h是步长。

以此类推，我们计算k3和k4分别在x0+h/2和x0+h处的斜率。

然后，我们通过加权组合这些斜率来计算下一个点的函数值y1：y1=y0+(h/6)*(k1+2*k2+2*k3+k4)。

3.算法步骤以下是使用四阶龙格库塔法求解常微分方程的详细步骤：步骤1：给定初始条件 y(x0) = y0，选择步长 h 和积分终点 x_end。

步骤2：计算积分步数n：n = (x_end - x0)/h。

步骤3：设置变量x=x0和y=y0，并将步长设置为h。

步骤4：对每个步数i=1,2,...,n，执行以下计算：4.1：计算斜率k1=f(x,y)。

4.2：计算斜率k2=f(x+h/2,y+h/2*k1)。

4.3：计算斜率k3=f(x+h/2,y+h/2*k2)。

4.4：计算斜率k4=f(x+h,y+h*k3)。

4.5：计算下一个点的函数值y1=y+(h/6)*(k1+2*k2+2*k3+k4)。

4.6：将x更新为x+h，y更新为y1步骤5：重复步骤4，直到达到积分终点 x_end。

4阶经典Runger-kutta格式解常微分方程

4阶经典Runger-kutta格式解常微分方程用Matlab软件:4阶经典Runger-kutta格式解常微分方程,y,f(x,y)考虑一阶常微分方程初值问题由Lagrange微分中值定理知道，存在y(x),y00,,(x,x)使得: nn，1yxyx,()(),n，1n,,y， ,,,()y(x),y(x)，h*y(,),y(x)，h*Kn1nn，h*,K,y(,),f(,,y(,))[x,x]其中称为在上的平均斜率 y(x)nn，1*[x,x]现在问题转化为如何对进行数值计算，可取在若干个点的斜率，或预Ky(x)nn，1rr,aKa,1报斜率值K,K,,,K的加权平均 ()作为的近似值。

设计K,,iii12ri,1,1i[x,x]K,K,,,K上若干个点斜率值或预报斜率值及加权系数使得差分格式 nn，112rry(x),y(x)，h*aK达到r阶，称为r阶Runge-kutta格式。

,，1nnii,1i 由此导出四阶经典Runge-kutta格式:y,y，h/6*(K，2K，2K，K)n，1n1234K,f(x，y)1nnK,f(x,y，h/2*K) 2n，1/2n1K,f(x,y，h/2*K)3n，1/2n2K,f(x,y，h/2*K)4n，1/2n3程序:(1)调用函数程序:function [x,y]=nark4(dyfun,xspan,y0,h) %用途:四阶经典Runge-Kutta格式解常微分方程y'=f(x,y),y(x0)=y0; %格式:[x,y]=nark4(dyfun,xspan,y0,h) %dyfun为函数f(x,y),xspan为求解区间[x0,xN],y0为初始值y(x0),h为步长，x 返回节点，y返回函数值解x=xspan(1):h:xspan(2);y(1)=y0;for n=1:length(x)-1k1=feval(dyfun,x(n),y(n));k2=feval(dyfun,x(n)+h/2,y(n)+h/2*k1);k3=feval(dyfun,x(n)+h/2,y(n)+h/2*k2);k4=feval(dyfun,x(n+1),y(n)+h*k3);y(n+1)=y(n)+h*(k1+2*k2+2*k3+k4)/6;endx=x';y=y';编写程序:>> dyfun=inline('y-2*x/y');>> [x,y]=nark4(dyfun,[0,1],1,0.1);plot(x,y,'-r') >> hold on>> grid on>> [x,y]=nark4(dyfun,[0,1],1,0.05); >> plot(x,y,'-g')>> hold on>> [x,y]=nark4(dyfun,[0,1],1,0.2);plot(x,y,'-b')。

微分方程组的龙格库塔公式求解matlab版

微分方程组的龙格-库塔公式求解matlab 版南京大学王寻1.一阶常微分方程组考虑方程组()()()()⎩⎨⎧====0000z x z ,z ,y ,x g 'z y x y ,z ,y ,x f 'y 其经典四阶龙格-库塔格式如下：对于n =0,1,2,...,计算()()⎪⎩⎪⎨⎧++++=++++=++4321143211226226L L L L h z z K K K K h y y n n n n 其中()()()()⎪⎪⎪⎩⎪⎪⎪⎨⎧+++=+++=⎪⎭⎫ ⎝⎛+++=⎪⎭⎫ ⎝⎛+++=⎪⎭⎫ ⎝⎛+++=⎪⎭⎫ ⎝⎛+++===33433422322311211211222222222222hL z ,hK y ,h x g L ,hL z ,hK y ,h x f K hL z ,hK y ,h x g L ,hL z ,hK y ,h x f K hL z ,hK y ,h x g L ,hL z ,hK y ,h x f K z ,y ,x g L ,z ,y ,x f K n n n n n n n n n n n n n n n n n n n n n n n n 下面给出经典四阶龙格-库塔格式解常微分方程组的matlab 通用程序：%marunge4s.m%用途：4阶经典龙格库塔格式解常微分方程组y'=f(x,y),y(x0)=y0%格式：[x,y]=marunge4s(dyfun,xspan,y0,h)%dyfun 为向量函数f(x,y),xspan 为求解区间[x0,xn],%y0为初值向量，h 为步长，x 返回节点，y 返回数值解向量function [x,y]=marunge4s(dyfun,xspan,y0,h)x=xspan(1):h:xspan(2);y=zeros(length(y0),length(x));y(:,1)=y0(:);for n=1:(length(x)-1)k1=feval(dyfun,x(n),y(:,n));k2=feval(dyfun,x(n)+h/2,y(:,n)+h/2*k1);k3=feval(dyfun,x(n)+h/2,y(:,n)+h/2*k2);k4=feval(dyfun,x(n+1),y(:,n)+h*k3);y(:,n+1)=y(:,n)+(h/6).*(k1+2*k2+3*k3+k4);end如下为例题：例1：取h =0.02，利用程序marunge4s.m 求刚性微分方程组()()⎩⎨⎧=-==--=10100209999010z ,z 'z ,y ,z .y .'y 的数值解，其解析解为：xx x .e z ,e e y 100100010---=+=解：首先编写M 函数dyfun.m%dyfun.m function f=dyfun(t,y)f(1)=-0.01*y(1)-99.99*y(2);f(2)=-100*y(2);f=f(:);然后编写一个执行函数：close all ;clear all ;clc;[x,y]=marunge4s(@dyfun,[0500],[21],0.002);plot(x,y);axis([-50500-0.52]);text(120,0.4,'y(x)');text(70,0.1,'z(x)');title('数值解');figure;y1=exp(-0.01*x)+exp(-100*x);%解析解，用来做对比的。

三阶、四阶龙格库塔函数matlab代码

三阶、四阶龙格库塔函数m a t l a b代码本页仅作为文档页封面，使用时可以删除This document is for reference only-rar21year.March三阶龙格—库塔法的计算公式为： )4(6)2,()2,2(),(3211213121K K K h y y hK hK y h x g K K h y h x g K y x g K i i i i i i i i +++=+-+=++==+ 三阶龙格—库塔公式的Matlab 程序代码：function y = DELGKT3_kuta(f, h,a,b,y0,varvec)format long;N = (b-a)/h;y = zeros(N+1,1);y(1) = y0;x = a:h:b;var = findsym(f);for i=2:N+1K1 = Funval(f,varvec,[x(i-1) y(i-1)]);K2 = Funval(f,varvec,[x(i-1)+h/2 y(i-1)+K1*h/2]);K3 = Funval(f,varvec,[x(i-1)+h y(i-1)-h*K1+K2*2*h]);y(i) = y(i-1)+h*(K1+4*K2+K3)/6;endformat short;DELGKT3_kuta函数运行时需要调用下列函数：function fv=Funval(f, varvec, varval)var= findsym(f);if length(var)<4if var(1)==varvec(1)fv=subs(f,varvec(1),varval(1));elsefv=subs(f,varvec(2),varval(2));endelsefv=subs(f,varvec,varval);end三阶龙格—库塔求解一阶常微分方程应用实例。

用三阶龙格—库塔法求下面常微分方程的数值解。

matlab用经典龙格库塔法求微分方程组初值问题程序

matlab用经典龙格库塔法求微分方程组初值问题程序经典龙格-库塔法是一种数值求解常微分方程的方法。

以下是一个使用MATLAB实现经典龙格-库塔法求解微分方程组的示例代码：```matlabfunction [t, y] = runge_kutta(f, y0, tspan, N)% f: 微分方程右边的函数句柄% y0: 初始值% tspan: 时间范围 [t0, tf]% N: 步数t = linspace(tspan(1), tspan(2), N+1); % 时间向量y = zeros(size(t)); % 初始化解向量y(1) = y0; % 设置初始值for i = 1:N% 计算四个点上的值k1 = f(t(i), y(i));k2 = f(t(i) + h/2, y(i) + h/2k1);k3 = f(t(i) + h/2, y(i) + h/2k2);k4 = f(t(i) + h, y(i) + hk3);% 更新解向量y(i+1) = y(i) + h/6(k1 + 2k2 + 2k3 + k4);endend```在上述代码中，我们定义了一个名为 `runge_kutta` 的函数，该函数接受微分方程右边的函数句柄 `f`、初始值 `y0`、时间范围 `tspan` 和步数 `N` 作为输入，并返回时间向量 `t` 和解向量 `y`。

在函数内部，我们首先生成时间向量 `t`，然后初始化解向量 `y`，并设置初始值 `y0`。

接下来，我们使用一个循环来迭代计算每个时间点上的值，并使用龙格-库塔公式更新解向量。

最后，我们返回时间向量 `t` 和解向量 `y`。

要使用该函数求解微分方程组，可以按照以下步骤进行：1. 定义微分方程右边的函数句柄 `f`，该函数接受时间和解向量作为输入，并返回微分方程的右侧值。

2. 定义初始值 `y0`。

3. 定义时间范围 `tspan`，该向量包含时间范围的起始和终止值。

(完整word版)龙格库塔法求微分方程matlab

(完整word 版)龙格库塔法求微分方程matlab龙格—库塔方法求解微分方程初值问题（数学1201+41262022+陈晓云）初值问题：y x x -+=2dxdy ，10≤≤x 1)0(y = 四阶龙格-库塔公式：()y x K n n ,f 1=⎪⎪⎪⎪⎭⎫ ⎝⎛+=+K h y x K n h n 122f ,2 ⎪⎭⎫ ⎝⎛++=K y x fK h n h n 232,2 ()K h y h x f K n n 34,++=()K K K K y y h n 43211n 226++++=+ 程序：1）建立四阶龙格-库塔函数function [ x,y ] = nark4( dyfun,xspan,y0,h )% dyfun 为一阶微分方程的函数；y0为初始条件；xspan 表示x 的区间；h 为区间的步长； x=xspan(1):h:xspan(2);y(1)=y0;for n=1:length(x)-1k1=feval(dyfun,x(n),y(n));k2=feval(dyfun,x(n)+h/2,y(n)+h/2*k1);(完整word版)龙格库塔法求微分方程matlab k3=feval(dyfun,x(n)+h/2,y(n)+h/2*k2);k4=feval(dyfun,x(n+1),y(n)+h*k3);y(n+1)=y(n)+h*(k1+k2*2+2*k3+k4)/6;endx=x;y=y;2）执行程序（m文件）dyfun=inline('x^2+x-y');[x,y1]=nark4(dyfun,[0,1],1,0.1);x=0:0.1:1;Format longy2=x.^2-x+1R4=y2-y1[x',y1',y2',R4']y2=dsolve('Dy=x^2+x-y','y(0)=1','x')plot(x,y1,'b*-')hold ony3=inline('x^2-x+1')fplot(y3,[0,1],'ro-')legend('R-K4','解析解')3）执行结果ans =X RK4近似值解析值0 1.000000000000000 1.000000000000000 0.100000000000000 0.910000208333333 0.910000000000000 0.200000000000000 0.840000396841146 0.840000000000000 0.300000000000000 0.790000567410084 0.790000000000000 0.400000000000000 0.760000721747255 0.760000000000000 0.500000000000000 0.750000861397315 0.750000000000000 0.600000000000000 0.760000987757926 0.760000000000000 0.700000000000000 0.790001102093746 0.790000000000000 0.800000000000000 0.840001205549083 0.8400000000000000.900000000000000 0.910001299159352 0.9100000000000001.000000000000000 1.000001383861433 1.000000000000000误差-0.000000208333333-0.000000396841146-0.000000567410084-0.000000721747255-0.000000861397315-0.000000987757926-0.000001102093746-0.000001205549083-0.000001299159352-0.000001383861433 y2 =x^2 - x + 1结果分析：初值问题的解析解为Y=x^2 - x + 1由图看出龙格库塔方法误差很小，具有很高的精度。

四阶龙格库塔法求解微分方程组

四阶龙格库塔法求解微分方程组四阶龙格-库塔法（RK4）是一种常用的数值求解微分方程组的方法。

它通过将微分方程转化为差分方程，通过逐步迭代来求得方程的近似解。

首先，考虑一个一阶微分方程组：dx/dt = f(x,t)其中，x是一个向量，f是一个向量函数，t是时间。

我们的目标是求解x(t)。

RK4法的基本思想是，在每个步骤中，通过预测下一个时间点的解并进行修正来逼近实际的解。

具体步骤如下：1.给定初始条件x(t0)和步长h，计算t=t0+h。

2. 计算k1 = hf(x(t0),t0)，即利用初始点(x(t0),t0)计算出的斜率。

3. 计算k2 = hf(x(t0)+k1/2,t0+h/2)，即利用中间点(x(t0)+k1/2,t0+h/2)计算出的斜率。

4. 计算k3 = hf(x(t0)+k2/2,t0+h/2)，同理得到另一个中间点的斜率。

5. 计算k4 = hf(x(t0)+k3,t0+h)，即利用最后一个中间点(x(t0)+k3,t0+h)计算的斜率。

6.根据k1，k2，k3和k4计算下一个点的值：x(t1)=x(t0)+(k1+2k2+2k3+k4)/67.重复步骤2到6，直到达到所需的终止时间。

接下来，让我们通过一个具体的例子来说明RK4方法的应用。

假设我们要求解如下的一阶微分方程组：dx/dt = f(x,y,z)dy/dt = g(x,y,z)dz/dt = h(x,y,z)其中f，g和h是已知的向量函数，我们需要找到x(t)，y(t)和z(t)。

首先，我们需要将该微分方程组转化为差分方程组。

我们将离散时间点t0，t1，t2，...，tn代入微分方程，得到：x(t0+h)=x(t0)+h*f(x(t0),y(t0),z(t0))y(t0+h)=y(t0)+h*g(x(t0),y(t0),z(t0))z(t0+h)=z(t0)+h*h(x(t0),y(t0),z(t0))然后，我们通过逐步迭代的方式来求解方程组。

matlab用四阶龙格库塔函数求解微分方程组

一、介绍Matlab作为一种强大的科学计算软件，提供了众多函数和工具来解决微分方程组。

其中，四阶龙格库塔函数是一种常用的数值方法，用于求解常微分方程组。

本文将介绍如何使用Matlab中的四阶龙格库塔函数来求解微分方程组，并对该方法的原理和实现进行详细说明。

二、四阶龙格库塔方法四阶龙格库塔方法是一种常用的数值方法，用于求解常微分方程组。

它是一种显式的Runge-Kutta方法，通过逐步逼近微分方程的解，在每一步使用多个中间值来计算下一步的解。

该方法通过四个中间值来计算下一步的状态，并且具有较高的精度和稳定性。

三、在Matlab中使用四阶龙格库塔方法求解微分方程组在Matlab中，可以使用ode45函数来调用四阶龙格库塔方法来解决微分方程组的问题。

ode45函数是Matlab提供的用于求解常微分方程组的函数，可以通过指定微分方程组以及初值条件来调用四阶龙格库塔方法来进行求解。

1. 定义微分方程组我们需要定义要求解的微分方程组。

可以使用Matlab中的匿名函数来定义微分方程组，例如：```matlabf = (t, y) [y(2); -sin(y(1))];```其中，f是一个匿名函数，用于表示微分方程组。

在这个例子中，微分方程组是y' = y2, y2' = -sin(y1)。

2. 指定初值条件和求解区间接下来，我们需要指定微分方程组的初值条件和求解区间。

初值条件可以通过指定一个初始时刻的状态向量来完成，例如：```matlabtspan = [0, 10];y0 = [0, 1];```其中，tspan表示求解区间，y0表示初值条件。

3. 调用ode45函数进行求解我们可以通过调用ode45函数来求解微分方程组的数值解。

具体的调用方式如下：```matlab[t, y] = ode45(f, tspan, y0);```其中，t和y分别表示求解的时间点和对应的状态值。

四、示例下面我们通过一个具体的例子来演示如何使用Matlab中的四阶龙格库塔方法来求解微分方程组。

matlab编的4阶龙格库塔法解微分方程的程序

matlab编的4阶龙格库塔法解微分方程的程序2010-03—10 20：16function varargout=saxplaxliu(varargin)clc，clearx0=0；xn=1。

2；y0=1；h=0。

1;[y,x］=lgkt4j（x0,xn，y0,h)；n=length（x)；fprintf（’ i x（i) y(i)\n’)；for i=1:nfprintf('%2d %4.4f ％4.4f\n’,i，x(i），y（i))；endfunction z=f（x，y）z=-2*x*y^2；function [y,x]=lgkt4j（x0，xn,y0，h)x=x0:h:xn；n=length（x）；y1=x;y1（1)=y0;for i=1:n-1K1=f(x(i)，y1(i））;K2=f（x（i)+h/2,y1（i)+h/2*K1)；K3=f(x(i)+h/2,y1（i)+h/2*K2);K4=f(x（i）+h,y1（i）+h*K3）;y1(i+1)=y1(i）+h/6＊(K1+2*K2+2＊K3+K4）;endy=y1；结果：i x（i） y（i)1 0.0000 1。

00002 0。

1000 0.99013 0。

2000 0.96154 0.3000 0。

91745 0.4000 0。

86216 0.5000 0。

80007 0。

6000 0。

73538 0。

7000 0.67119 0。

8000 0。

609810 0。

9000 0。

552511 1.0000 0.500012 1.1000 0.452513 1.2000 0.4098yi+1=yi+h＊（ K1+ K2)/2K1=f（xi,yi)K2=f(xi+h，yi+h*K1)依次类推，如果在区间［xi，xi+1］内多预估几个点上的斜率值K1、K2、……Km,并用他们的加权平均数作为平均斜率K*的近似值，显然能构造出具有很高精度的高阶计算公式。