基于matlab的龙格库塔法求解布拉修斯方程
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是初始条件。
matlab-欧拉方法和龙格库塔方法的小实例

题一:a)y’=y+2x , 欧拉方法:112()2n n h y y k k +=++,12n n k y x =+,2112()n n k y hk x +=++; 龙格-库塔方法:11234(22)6n n h y y k k k k +=++++,12n n k y x =+,12222n n k h k y h x ⎛⎫=+++ ⎪⎝⎭,23222n n k h k y h x ⎛⎫=+++ ⎪⎝⎭,432()n n k y hk x h =+++ 精确解:y=3e x -2x-2。
以步长h=0.1 在0<=x<=1内的计算结果如下所示:0.1000 1.1150 1.1155 1.11550.2000 1.2631 1.2642 1.26420.3000 1.4477 1.4496 1.44960.4000 1.6727 1.6755 1.67550.5000 1.9423 1.9462 1.94620.6000 2.2613 2.2664 2.26640.7000 2.6347 2.6413 2.64130.8000 3.0684 3.0766 3.07660.9000 3.5685 3.5788 3.57881.0000 4.1422 4.1548 4.1548b)文案 编辑词条B 添加义项 ?文案,原指放书的桌子,后来指在桌子上写字的人。
现在指的是公司或企业中从事文字工作的职位,就是以文字来表现已经制定的创意策略。
文案它不同于设计师用画面或其他手段的表现手法,它是一个与广告创意先后相继的表现的过程、发展的过程、深化的过程,多存在于广告公司,企业宣传,新闻策划等。
基本信息中文名称文案外文名称Copy目录1发展历程2主要工作3分类构成4基本要求5工作范围6文案写法7实际应用折叠编辑本段发展历程汉字"文案"(wén àn)是指古代官衙中掌管档案、负责起草文书的幕友,亦指官署中的公文、书信等;在现代,文案的称呼主要用在商业领域,其意义与中国古代所说的文案是有区别的。
龙格库塔方法及其matlab实现

龙格-库塔方法及其matlab实现摘要:本文的目的数值求解微分方程精确解,通过龙格-库塔法,加以利用matlab为工具达到求解目的。
龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法,用于数值求解微分方程。
MatLab软件是由美国Mathworks公司推出的用于数值计算和图形处理的科学计算系统环境。
MatLab是英文MATrix LABoratory(矩阵实验室)的缩写。
在MratLab环境下,用户可以集成地进行程序设计、数值计算、图形绘制、输入输出、文件管理等各项操作。
关键词:龙格-库塔 matlab 微分方程1.前言1.1:知识背景龙格-库塔法(Runge-Kutta)是用于非线性常微分方程的解的重要的一类隐式或显式迭代法。
这些技术由数学家卡尔·龙格和马丁·威尔海姆·库塔于1900年左右发明。
通常所说的龙格库塔方法是相对四阶龙格库塔而言的,成为经典四阶龙格库塔法。
该方法具有精度高,收敛,稳定,计算过程中可以改变步长不需要计算高阶导数等优点,但是仍需计算在一些点上的值,比如四阶龙格-库塔法没计算一步需要计算四步,在实际运用中是有一定复杂性的。
Matlab是在20世纪七十年代后期的事:时任美国新墨西哥大学计算机科学系主任的Cleve Moler教授出于减轻学生编程负担的动机,为学生设计了一组调用LINPACK和EISPACK库程序的“通俗易用”的接口,此即用FORTRAN编写的萌芽状态的MATLAB。
经几年的校际流传,在Little的推动下,由Little、Moler、Steve Bangert合作,于1984年成立了MathWorks公司,并把MATLAB正式推向市场。
从这时起,MATLAB的内核采用C语言编写,而且除原有的数值计算能力外,还新增了数据图视功能。
MATLAB以商品形式出现后,仅短短几年,就以其良好的开放性和运行的可靠性,使原先控制领域里的封闭式软件包(如英国的UMIST,瑞典的LUND和SIMNON,德国的KEDDC)纷纷淘汰,而改以MATLAB为平台加以重建。
Matlab中龙格-库塔(Runge-Kutta)方法原理及实现

函数功能编辑本段回目录ode是专门用于解微分方程的功能函数,他有ode23,ode45,ode23s等等,采用的是Runge-Kutta算法。
ode45表示采用四阶,五阶runge-kutta单步算法,截断误差为(Δx)³。
解决的是Nonstiff(非刚性)的常微分方程.是解决数值解问题的首选方法,若长时间没结果,应该就是刚性的,换用ode23来解.使用方法编辑本段回目录[T,Y] = ode45(odefun,tspan,y0)odefun 是函数句柄,可以是函数文件名,匿名函数句柄或内联函数名tspan 是区间[t0 tf] 或者一系列散点[t0,t1,...,tf]y0 是初始值向量T 返回列向量的时间点Y 返回对应T的求解列向量[T,Y] = ode45(odefun,tspan,y0,options)options 是求解参数设置,可以用odeset在计算前设定误差,输出参数,事件等[T,Y,TE,YE,IE] =ode45(odefun,tspan,y0,options)在设置了事件参数后的对应输出TE 事件发生时间YE 事件解决时间IE 事件消失时间sol =ode45(odefun,[t0 tf],y0...)sol 结构体输出结果应用举例编辑本段回目录1 求解一阶常微分方程程序:一阶常微分方程odefun=@(t,y) (y+3*t)/t^2; %定义函数tspan=[1 4]; %求解区间y0=-2; %初值[t,y]=ode45(odefun,tspan,y0);plot(t,y) %作图title('t^2y''=y+3t,y(1)=-2,1<t<4')legend('t^2y''=y+3t')xlabel('t')ylabel('y')% 精确解% dsolve('t^2*Dy=y+3*t','y(1)=-2')% ans =一阶求解结果图% (3*Ei(1) - 2*exp(1))/exp(1/t) - (3*Ei(1/t))/exp(1/t)2 求解高阶常微分方程关键是将高阶转为一阶,odefun的书写.F(y,y',y''...y(n-1),t)=0用变量替换,y1=y,y2=y'...注意odefun方程定义为列向量dxdy=[y(1),y(2)....]程序:function Testode45tspan=[3.9 4.0]; %求解区间y0=[2 8]; %初值[t,x]=ode45(@odefun,tspan,y0);plot(t,x(:,1),'-o',t,x(:,2),'-*')legend('y1','y2')title('y'' ''=-t*y + e^t*y'' +3sin2t')xlabel('t')ylabel('y')function y=odefun(t,x)y=zeros(2,1); % 列向量y(1)=x(2);y(2)=-t*x(1)+exp(t)*x(2)+3*sin(2*t);endend高阶求解结果图相关函数编辑本段回目录ode23, ode45, ode113, ode15s, ode23s, ode23t, ode23tbMatlab中龙格-库塔(Runge-Kutta)方法原理及实现(自己写的,非直接调用)龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法。
4阶经典Runger-kutta格式解常微分方程

用Matlab 软件:4阶经典Runger-kutta 格式解常微分方程考虑一阶常微分方程初值问题00)(),(y x y y x f y =='由Lagrange 微分中值定理知道,存在),(1+∈n n x x ε使得:⇒-='+hx y x y y n n )()()(1ε, *++='+=K h x y y h x y x y n n n *)()(*)()(1ε 其中))(,()(*εεεy f y K ='=称为)(x y 在],[1+n n x x 上的平均斜率现在问题转化为如何对*K 进行数值计算,可取)(x y 在],[1+n n x x 若干个点的斜率,或预报斜率值r K K K ⋅⋅⋅21,的加权平均∑=r i i i K a 1 (∑==ri i a 11)作为*K 的近似值。
设计],[1+n n x x 上若干个点斜率值或预报斜率值 及加权系数r K K K ⋅⋅⋅21,使得差分格式 ∑=++=ri i i n n K a h x y x y 11*)()(达到r 阶,称为r 阶Runge-kutta 格式。
由此导出四阶经典Runge-kutta 格式:)*2/,()*2/,()*2/,()()22(*6/32/1422/1312/12143211K h y x f K K h y x f K K h y x f K y x f K K K K K h y y n n n n n n n n n n +=+=+=+=++++=++++程序:(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.龙格库塔法的基本原理2.龙格库塔法的发展历程三、MATLAB 实现四阶龙格库塔法1.MATLAB 中龙格库塔法的函数2.四阶龙格库塔法的MATLAB 实现四、龙格库塔法解方程组的应用1.线性方程组的求解2.非线性方程组的求解五、结论正文:一、引言在数学领域,求解方程组是一项基本任务。
龙格库塔法作为高效数值求解线性方程组的方法,被广泛应用于各个领域。
MATLAB 作为一款强大的数学软件,可以方便地实现龙格库塔法求解方程组。
本文将介绍MATLAB 中四阶龙格库塔法解方程组的原理、实现与应用。
二、龙格库塔法介绍1.龙格库塔法的基本原理龙格库塔法(Runge-Kutta method)是一种求解常微分方程初值问题的数值方法。
它通过求解一组线性方程来逼近微分方程的解,具有较高的数值稳定性和精度。
龙格库塔法可以分为四阶、五阶等多种形式,其中四阶龙格库塔法是较为常用的一种。
2.龙格库塔法的发展历程龙格库塔法由德国数学家卡尔·龙格(Carl Runge)和英国数学家詹姆斯·库塔(James Kutta)分别在1900 年和1901 年独立发现。
自那时以来,龙格库塔法在数学、物理、工程等领域得到了广泛应用,并发展出了多种改进和扩展。
三、MATLAB 实现四阶龙格库塔法1.MATLAB 中龙格库塔法的函数在MATLAB 中,可以使用内置函数ode45、ode23、ode113 等实现龙格库塔法求解常微分方程。
这些函数分别对应四阶、五阶和三阶龙格库塔法。
2.四阶龙格库塔法的MATLAB 实现以下是一个使用MATLAB 实现四阶龙格库塔法求解方程组的示例:```matlabfunction [x, status] = solve_system_with_ode45(A, B, x0)% 定义方程组func = @(t, x) A * x + B;% 初始条件x0 = [1; 2];% 时间区间tspan = [0, 10];% 求解[x, status] = ode45(func, tspan, x0);end```四、龙格库塔法解方程组的应用1.线性方程组的求解线性方程组在数学、物理、工程等领域具有广泛应用。
matlab编写龙格库塔法或欧拉法求解常微分方程数值解

龙格库塔法(Runge-Kutta method)和欧拉法(Euler's method)是两种常用的数值求解常微分方程的方法。
这里分别给出它们的MATLAB实现:1. 龙格库塔法(Runge-Kutta method):```matlabfunction [y, t] = runge_kutta(f, y0, t0, tf, h)% f: 微分方程函数,输入为[y, t],输出为dy/dt% y0: 初始值% t0: 初始时间% tf: 结束时间% h: 步长N = round((tf - t0) / h); % 计算迭代次数t = zeros(1, N + 1); % 初始化时间向量y = zeros(size(y0), N + 1); % 初始化解向量t(1) = t0;y(:, 1) = y0;for i = 1:Nk1 = h * f(y(:, i), t(i));k2 = h * f(y(:, i) + k1 / 2, t(i) + h / 2);k3 = h * f(y(:, i) + k2 / 2, t(i) + h / 2);k4 = h * f(y(:, i) + k3, t(i) + h);y(:, i + 1) = y(:, i) + (k1 + 2 * k2 + 2 * k3 + k4) / 6;t(i + 1) = t(i) + h;endend```2. 欧拉法(Euler's method):```matlabfunction [y, t] = euler_method(f, y0, t0, tf, h)% f: 微分方程函数,输入为[y, t],输出为dy/dt% y0: 初始值% t0: 初始时间% tf: 结束时间% h: 步长N = round((tf - t0) / h); % 计算迭代次数t = zeros(1, N + 1); % 初始化时间向量y = zeros(size(y0), N + 1); % 初始化解向量t(1) = t0;y(:, 1) = y0;for i = 1:Ny(:, i + 1) = y(:, i) + h * f(y(:, i), t(i));t(i + 1) = t(i) + h;endend```使用这两个函数时,需要定义一个表示微分方程的函数`f`,例如:```matlabfunction dydt = my_ode(y, t)dydt = -y; % 一个简单的一阶线性微分方程:dy/dt = -yend```然后调用相应的求解函数,例如:```matlaby0 = 1; % 初始值t0 = 0; % 初始时间tf = 5; % 结束时间h = 0.1; % 步长[y_rk, t_rk] = runge_kutta(@my_ode, y0, t0, tf, h);[y_euler, t_euler] = euler_method(@my_ode, y0, t0, tf, h);```。
基于matlab的龙格库塔法求解布拉修斯方程

Runge —Kutta 法求解布拉修斯解摘要薄剪切层方程主要有三种解法,即相似解,非相似条件下对偏微分方程组的数值解和近似解。
布拉修斯解是布拉修斯于1908年求出的,它是零攻角沿平板流动的相似解。
本文用四阶Runge —Kutta 法求解高阶微分方程的方法,并用matlab 编程实现,求得了与实际层流边界层相符合的数值解。
关键词:布拉修斯解,相似解,Runge —Kutta 法,数值解。
1 布拉修斯近似解方程二维定常不可压缩层流边界层的方程为:0=∂∂+∂∂yvx u (1)22yuv dx d y u v x u u u u e e ∂∂+=∂∂+∂∂(2)边界条件为:0=y )(,0x v u v w ==:δ=y )(x u u e =将式(1)和式(2)进行法沃克纳—斯坎变换(简称F —S 变换),将边界层方程无量纲化,即设(3)y x v u e 5.0⎪⎪⎭⎫⎝⎛=η(4)x x =得出F —S 变换后的动量方程(5)()[]()[]⎪⎭⎫ ⎝⎛∂∂''-∂'∂'='-+''++'''+x f f xf f x f m f f m f t k221211其中k 为流动类型指标,横曲率项t 为(6)212120cos 211⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡⎪⎪⎭⎫ ⎝⎛⎪⎪⎭⎫ ⎝⎛++-=ηφe u vxL r L t m 是量纲一的压力梯度参数,定义为(7)xd du u x m ee =其边界条件变为:0=η0='f:∞=η1='f 对于二维平面实壁流动()可以忽略横曲率项t 的轴对称流动,式(5)成:0=η0=w f 为(8)()[]⎪⎭⎫ ⎝⎛∂∂''-∂'∂'='-+''++'''x f f xf f x f m f f m f 2121根据相似解的定义,方程(8)中的函数f 若式相似的,则它应只与η有关而与x 无关,即对x 的偏导数应为零。
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
Runge —Kutta 法求解布拉修斯解摘要薄剪切层方程主要有三种解法,即相似解,非相似条件下对偏微分方程组的数值解和近似解。
布拉修斯解是布拉修斯于1908年求出的,它是零攻角沿平板流动的相似解。
本文用四阶Runge —Kutta 法求解高阶微分方程的方法,并用matlab 编程实现,求得了与实际层流边界层相符合的数值解。
关键词:布拉修斯解,相似解,Runge —Kutta 法,数值解。
1 布拉修斯近似解方程二维定常不可压缩层流边界层的方程为:0=∂∂+∂∂yvx u (1)22yuv dx d y u v x u u u u e e ∂∂+=∂∂+∂∂ (2)边界条件为:0=y )(,0x v u vw==:δ=y )(x u u e =将式(1)和式(2)进行法沃克纳—斯坎变换(简称F —S 变换),将边界层方程无量纲化,即设y x v u e 5.0⎪⎪⎭⎫⎝⎛=η (3)x x = (4)得出F —S 变换后的动量方程()[]()[]⎪⎭⎫ ⎝⎛∂∂''-∂'∂'='-+''++'''+x f f xf f x f m f f m f t k221211 (5)其中k 为流动类型指标,横曲率项t 为212120cos 211⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡⎪⎪⎭⎫ ⎝⎛⎪⎪⎭⎫ ⎝⎛++-=ηφe u vx L r L t (6) m 是量纲一的压力梯度参数,定义为xd du u x m ee =(7)其边界条件变为:0=η 0='f:∞=η 1='f对于二维平面实壁流动(:0=η0=w f )可以忽略横曲率项t 的轴对称流动,式(5)成为()[]⎪⎭⎫ ⎝⎛∂∂''-∂'∂'='-+''++'''x f f xf f x f m f f m f 2121 (8) 根据相似解的定义,方程(8)中的函数f 若式相似的,则它应只与η有关而与x 无关,即对x 的偏导数应为零。
于是方程(8)应成为()[]01212='-+''++'''f m f f m f (9) 若f w 为常数,则方程(9)的边界条件为:0=η 常数==w f f ; 0='='wf f :∞=η 1='f2 布拉修斯解布拉修斯于1908年求出了零攻角沿平板流动的解。
这时 0==m u e 常数; 因而方程(9)成为021=''+'''f f f (10) 此即布拉修斯方程。
对于实壁,0=w f ,边界条件成为:0=η 0==w f f ; 0='='wf f :∞=η 1='f3 Runge —Kutta 法求解Runge —Kutta 通过将高阶微分方程化为一阶线性方程组,从而解出高阶方程的数值解。
在方程(10)中令η=0f⎪⎪⎪⎩⎪⎪⎪⎨⎧=====ηηηηd df d f d f d df d df f f f 2221213 (11) 于是方程(10)变为⎪⎪⎪⎩⎪⎪⎪⎨⎧-======313210333321022232101121),,,(),,,(),,,(f f f f f f T d df f f f f f T d df f f f f f T d df ηηη (12) 当区步长为h ,有四阶Runge —Kutta 的形式如下()⎪⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎨⎧++++=++++=++++==++++=+),,,()2,2,2,2()2,2,2,2(),,,(2263,3,33,2,23,1,13,0,04,2,3,32,2,22,1,12,0,03,1,3,31,2,21,1,11,0,02,,3,2,1,01,4,3,2,1,,1,hK f hK f hK f hK f T K K h f K h f K h f K h f T K K h f K h f K h f K h f T K f f f f T K K K K K h f f i i i i i i i i i i i i i i i i i i i i i i i i i i i i j i j i (13)使用matlab 软件取步长为0.2,迭代100步视作η→无穷大。
迭代到第40步左右就收敛了,迭代结果如下(本文附录有全程序源代码)表格 1平板层流边界层方程的数值解1.20.237950.393780.316591.40.322990.456270.307871.60.420330.516760.296671.80.529520.574760.2829320.650030.629770.266752.20.78120.681320.248352.40.92230.728990.228092.6 1.07250.772460.206462.8 1.2310.811510.184013 1.39680.846050.161363.2 1.56910.876090.139133.4 1.7470.901770.117883.6 1.92950.923330.0980873.8 2.1160.941120.0801264 2.30580.955520.0642354.2 2.49810.966960.0505214.4 2.69240.975880.0389744.6 2.88830.982690.0294854.8 3.08530.987790.0218735 3.28330.991550.0159085.2 3.48190.994250.0113435.4 3.68090.996160.0079295.6 3.88030.997480.0054345.8 4.07990.998380.003656 4.27960.998980.0024036.2 4.47950.999370.0015516.4 4.67940.999620.0009826.6 4.87930.999770.0006096.8 5.07930.999870.000377 5.27930.999930.0002217.2 5.47930.999960.0001297.4 5.67930.999987.38E-057.6 5.87930.99999 4.15E-057.8 6.07931 2.28E-058 6.27931 1.23E-058.2 6.47931 6.52E-068.4 6.67931 3.38E-068.6 6.87931 1.72E-068.87.079318.59E-0797.27931 4.20E-07由上表可以看出,数值解与布拉修斯解符合程度相当好。
4 结语(1) 布拉修斯用相似变换将N—S方程简化,简化为一维微分方程求解并得到了与实际层流现象相符合的结果。
(2)Runge—Kutta方法用来求解高阶微风方程,有相当高的精度。
参考文献:[1] 陈懋章.粘性流体力学基础.北京.高等教育出版社,2002.[2] 复旦大学数学系《微分方程及其数值解》编写组.《微分方程及其数值解》.上海.复旦大学出版社.2001年[3] 清华大学工程力学系编.流体力学基础:上册,下册.北京:机械工业出版社,1980年附录源程序代码如下% this file is to settle the bulaxiusi function with the method of% Runge-Kutta.% f1 stands for the function f% f2 stands for the function df1/du% f3 stands for the function df2/duf(1:3,1:100)=0; %取A的1,3行,1,100列。
f(3,1)=0.33206;x(1:101)=0;% h stands for the step length;h=0.2;% k1, k2,k3,k4 stands for the coefficient of Rung_Kutta of f1k=[0,0,0,0;0,0,0,0;0,0,0,0];for i=1:100;k(1,1)=f(2,i);k(2,1)=f(3,i);k(3,1)=-f(1,i)*f(3,i)/2k(1,2)=f(2,i)+h*k(2,1)/2;k(2,2)=f(3,i)+h*k(3,1)/2;k(3,2)=-(f(1,i)+h*k(1,1)/2)*(f(3,i)+h*k(3,1)/2)/2;k(1,3)=f(2,i)+h*k(2,2)/2;k(2,3)=f(3,i)+h*k(3,2)/2;k(3,3)=-(f(1,i)+h*k(1,2)/2)*(f(3,i)+h*k(3,2)/2)/2;k(1,4)=f(2,i)+h*k(2,3);k(2,4)=f(3,i)+h*k(3,3);k(3,4)=-(f(1,i)+h*k(1,3))*(f(3,i)+h*k(3,3))/2;f(1,i+1)=f(1,i)+h*(k(1,1)+2*k(1,2)+2*k(1,3)+k(1,4))/6;f(2,i+1)=f(2,i)+h*(k(2,1)+2*k(2,2)+2*k(2,3)+k(2,4))/6;f(3,i+1)=f(3,i)+h*(k(3,1)+2*k(3,2)+2*k(3,3)+k(3,4))/6;x(i+1)=x(i)+h;end。