四阶龙格库塔法解微分方程教学内容
- 格式:docx
- 大小:98.25 KB
- 文档页数:7
四阶龙格库塔法解微分方程一、四阶龙格库塔法解一阶微分方程①一阶微分方程:cos y t ,初始值y(0)=0,求解区间[0 10]。
MATLAB 程序:%%%%%%%%%%% 四阶龙哥库塔法解一阶微分方程%%%%%%%%%%% y'=cost%%%%%%%%%%% y(0)=0, 0≤t ≤10,h=0.01%%%%%%%%%%% y=sinth=0.01;hf=10;t=0:h:hf;y=zeros(1,length(t));y(1)=0;F=@(t,y)(cos(t));for i=1:(length(t)-1)k1=F(t(i),y(i));k2=F(t(i)+h/2,y(i)+k1*h/2);k3=F(t(i)+h/2,y(i)+k2*h/2);k4=F(t(i)+h,y(i)+k3*h);y(i+1)=y(i)+1/6*(k1+2*k2+2*k3+k4)*h;endsubplot(211)plot(t,y,'-')xlabel('t');ylabel('y');title('Approximation');span=[0,10];p=y(1);[t1,y1]=ode45(F,span,p);subplot(212)plot(t1,y1)xlabel('t');ylabel('y');title('Exact');图1②一阶微分方程:()22*/x t x x t =- ,初始值x(1)=2,求解区间[1 3]。
MATLAB 程序: %%%%%%%%%%% 四阶龙哥库塔法解微分方程%%%%%%%%%%% x'(t)=(t*x-x^2)/t^2%%%%%%%%%%% x(1)=2, 1≤t ≤3, h=1/128%%%%%%%%%%% 精确解:x(t)=t/(0.5+lnt)h=1/128; %%%%% 步长tf=3;t=1:h:tf;x=zeros(1,length(t));x(1)=2; %%%%% 初始值F_tx=@(t,x)(t.*x-x.^2)./t.^2;for i=1:(length(t)-1)k_1=F_tx(t(i),x(i));k_2=F_tx(t(i)+0.5*h,x(i)+0.5*h*k_1);k_3=F_tx((t(i)+0.5*h),(x(i)+0.5*h*k_2));k_4=F_tx((t(i)+h),(x(i)+k_3*h));x(i+1)=x(i)+(1/6)*(k_1+2*k_2+2*k_3+k_4)*h; endsubplot(211)plot(t,x,'-');xlabel('t');ylabel('x');legend('Approximation');%%%%%%%%%%%%%%%%%%%%%%%%%%%% ode45求精确解t0=t(1);x0=x(1);xspan=[t0 tf];[x_ode45,y_ode45]=ode45(F_tx,xspan,x0);subplot(212)plot(x_ode45,y_ode45,'--');xlabel('t');ylabel('x');legend('Exact');图2二、四阶龙格库塔法解二阶微分方程①二阶微分方程:cos y t ,初始值y(0)=0,y'(0)=-1,求解区间[0 10]。
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。
四阶龙格库塔法(Runge-Kutta )求解微分方程张晓颖(天津大学 材料学院 学号:1012208027)1 引言计算传热学中通常需要求解常微分方程。
这类问题的简单形式如下:{),(')(00y x f y y x y == (1)虽然求解常微分方程有各种各样的解析方法,但解析方法只能用来求解一些特殊类型的方程,实际问题中的多数微分方程需要采用数值解法求解。
初值问题(1)的数值解法有个基本特点,它们采取“步进式”,即求解过程顺着节点排序一步一步向前推进。
这类算法是要给出用已知信息y n 、 y n −i ……计算y n +1的递推公式即可。
2 龙格库塔法(Runge-Kutta )介绍假设对于初值问题(1)有解 y = y (x ) ,用 Taylor 展开有:......)(!3)(!2)()()(321+'''+''+'+=+n n n n n x y h x y h x y h x y x y (2) 龙格库塔法(Runge-Kutta )实质上是间接的使用 Taylor 级数法的一种方法。
对于差商hx y x y n n )()(1-+,根据微分中值定理,存在 0 < θ < 1 ,使得:)()()(1h x y hx y x y n n θ+'=-+ (3)于是对于 y = y (x ) ,可得到:))(,()()(1h x y h x hf x y x y n n n n θθ+++=+ (4)设))(,(*h x y h x f K n n θθ++=,为区间 [x n , x n +1 ] 上的平均斜率。
四阶龙格库塔格式中的*K 由下式计算得到:⎪⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎨⎧++=++=++==++++=+),()2,2()2,2(),()22(6342312143211hK y h x f K K h y h x f K K h y h x f K y x f K K K K K h y y n n n n nn n n n n (5) 四阶龙格库塔法(Runge-Kutta )的每一步需要四次计算函数值f ,其截断误差更低,计算的精度更高。
四阶龙格库塔解二阶微分方程四阶龙格库塔法是一种常用的数值解微分方程的算法,适用于解一阶、二阶以及高阶微分方程。
本文将围绕“四阶龙格库塔解二阶微分方程”展开,并分步骤进行阐述。
第一步:转换为一阶微分方程组对于二阶微分方程 y''=f(x,y,y'),可以将其转换为一阶微分方程组。
令 u=y',则 y''=du/dx,原方程变为 du/dx=f(x,y,u),整理得到以下一阶微分方程组:dy/dx=udu/dx=f(x,y,u)第二步:确定步长和初始条件在使用四阶龙格库塔法时,需要确定一个步长 h,以及初始条件y0 和 u0。
步长 h 决定了每一步积分时横坐标 x 的变化量。
初始条件 y0 和 u0 分别代表了 x=0 时的纵坐标和导数值。
第三步:进行数值积分按照四阶龙格库塔法的步骤进行数值积分,可得到第 k 步的纵坐标 yk 和导数值 uk:k1 = hf(xk,yk,uk)j1 = hukk2 = hf(xk+h/2,yk+j1/2,uk+k1/2)j2 = h(uk+k1/2)k3 = hf(xk+h/2,yk+j2/2,uk+k2/2)j3 = h(uk+k2/2)k4 = hf(xk+h,yk+j3,uk+k3)j4 = h(uk+k3)yk+1 = yk+(j1+2j2+2j3+j4)/6uk+1 = uk+(k1+2k2+2k3+k4)/6其中,k1 到 k4 和 j1 到 j4 分别为由函数 f(x,y,u) 求得的中间变量。
yk 和 uk 即为求得的第 k 步的纵坐标和导数值,yk+1 和uk+1 分别为求得的第 k+1 步的纵坐标和导数值。
第四步:循环迭代将上述计算步骤循环执行至所需的精度或区间内积分完成为止。
即对于每一步 k,根据当前的 yk 和 uk 求解下一步的 yk+1 和 uk+1,并记录下这些变量的值。
至此,关于“四阶龙格库塔解二阶微分方程”的阐述就结束了。
中南大学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 圆。
计算机仿真作业2(院、系) 专业 班一、实验目的熟悉matlab 环境和基本操作。
二、实验内容熟悉matlab 环境及在自动控制中的应用。
面向微分方程的数字仿真1. 分别使用(1)四阶龙格-库塔解微分方程方法、(2)用matab 的ode45函数求解具有如下闭环传递函数的系统的阶跃响应。
43210()8364010s s s s s φ=++++解:(1)四阶龙格-库塔解微分方程方法r=1numo=[10];demo=[1 8 36 40 10];numh=0.002;demh=1;[num,den]=feedback(numo,demo,numh,demh);[A,B,C,D]=tf2ss(num,den);X=[zeros(length(A),1)];Y=0;t=0;Tf=5;h=0.02;n=Tf/h;for i=1:nK1=A*X+B*r;K2=A*(X+h*K1/2)+B*r;K3= A*(X+h*K2/2)+ B*r;K4=A*(X+h*K3)+B*r;X=X+h*(K1+2*K2+2*K3+K4)/6;Y=[Y,C*X+D*r];t=[t,t(i)+h];endplot(t,Y)(2) ode45函数解微分方程法:(2) ode45函数解微分方程法1、把连续函数转换成一阶微分方程组:用MATLAB编程>> num=[10];den=[1 8 36 40 10];[A B C D]=tf2ss(num,den)A =-8 -36 -40 -101 0 0 00 1 0 00 0 1 0B =1C =0 0 0 10D =0 所以模型为:function fun=fun(t,x)u=100;fun=[-8*x(1)-36*x(2)-40*x(3)-10*x(4)+u;x(1);x(2);x(3)]x0=[0,0,0,0],tspa=[0,5];x0 =0 0 0 0>> [t,y]=ode45('fun',tspa,x0)plot(t,y)[]x y u x x 10 0 0 000010 1 0 0 0 0 1 0 0 0 0 1 10- 40- 36- 8-=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡+⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=。
四阶龙格库塔方法求解n自由度二阶微分方程在数值计算中,龙格-库塔方法(Runge-Kutta method)是一种常用的求解常微分方程(ODE)的数值方法。
它是一种多步法(multiple-step method),通过使用不同的近似值来迭代计算ODE 的解。
其中,四阶龙格-库塔方法是最常用的龙格-库塔方法之一,也是一种高阶方法,可用于求解n自由度二阶微分方程。
下面将介绍四阶龙格-库塔方法的原理、步骤和应用。
首先,我们来回顾一下二阶微分方程的一般形式:y''(y)=y(y,y(y),y'(y))这里,y(y,y(y),y'(y))是已知的函数,y(y)是我们要求解的未知函数。
为了求解这个二阶微分方程,我们需要将其转化为一个一阶微分方程的求解问题。
令y(y)=y(y)y(y)=y'(y)这样,我们可以将原始的二阶微分方程转化为如下的一阶微分方程组:y'(y)=y(y)y'(y)=y(y,y(y),y(y))现在,我们可以利用四阶龙格-库塔方法来求解这个一阶微分方程组。
四阶龙格-库塔方法基于泰勒展开的思想,通过使用一系列的近似值来逼近方程的解。
该方法的一般形式为:y1=ℎy(yy,yy)y2=ℎy(yy+ℎ/2,yy+y1/2)y3=ℎy(yy+ℎ/2,yy+y2/2)y4=ℎy(yy+ℎ,yy+y3)yy+1=yy+1/6(y1+2y2+2y3+y4)yy+1=yy+ℎ其中,yy是当前时间步,yy+1是下一个时间步,yy是当前位置的近似解,yy+1是下一个位置的近似解,ℎ是时间步长,y1、y2、y3和y4是中间的近似值。
四阶龙格-库塔方法的步骤如下:Step 1:给定初始条件我们需要提供初始条件,包括初始位置y0和初始速度y0,以及时间步长ℎ、求解的时间区间[y0,yy]和离散时间步数y。
y0=y(y0)y0=y'(y0)Step 2:进行迭代计算从n=0开始,进行y步的迭代计算,其中y=(yy−y0)/ℎ。
四阶龙格库塔法求解微分方程组微分方程是数学中的一个重要分支,它描述了自然界中许多现象的发展规律。
在实际应用中,我们经常需要求解微分方程组,以得到系统的演化规律和性质。
本文将介绍一种常用的求解微分方程组的数值方法——四阶龙格库塔法。
一、微分方程组的数值解法微分方程组是形如下式的方程集合:begin{cases} frac{dy_1}{dx}=f_1(x,y_1,y_2,cdots,y_n) frac{dy_2}{dx}=f_2(x,y_1,y_2,cdots,y_n) cdotsfrac{dy_n}{dx}=f_n(x,y_1,y_2,cdots,y_n) end{cases} 其中,$y_1,y_2,cdots,y_n$是关于自变量$x$的未知函数,$f_1,f_2,cdots,f_n$是已知的函数。
求解微分方程组的数值方法主要有以下两种:1.欧拉法欧拉法是最简单的数值方法之一,它的基本思想是将微分方程组中的导数用差分代替,从而得到一个递推公式。
具体而言,欧拉法的递推公式为:$$y_{i+1}=y_i+hf(x_i,y_i)$$其中,$h$是步长,$x_i$和$y_i$分别是第$i$个点的自变量和因变量的值。
欧拉法的精度较低,容易出现数值误差,但是它的计算速度快,适用于一些简单的微分方程组。
2.龙格库塔法龙格库塔法是一种常用的高精度数值方法,它的基本思想是将微分方程组中的导数用一系列加权的差分代替,从而得到一个递推公式。
其中,四阶龙格库塔法是最为常用的一种。
具体而言,四阶龙格库塔法的递推公式为:begin{aligned} k_1&=hf(x_i,y_i)k_2&=hf(x_i+frac{h}{2},y_i+frac{k_1}{2})k_3&=hf(x_i+frac{h}{2},y_i+frac{k_2}{2})k_4&=hf(x_i+h,y_i+k_3)y_{i+1}&=y_i+frac{k_1+2k_2+2k_3+k_4}{6} end{aligned} 其中,$k_1,k_2,k_3,k_4$是加权的差分,$h$是步长,$x_i$和$y_i$分别是第$i$个点的自变量和因变量的值。
2013-2014(1)专业课程实践论文题目:四阶龙格—库塔法一、算法理论由定义可知,一种数值方法的精度与局部截断误差()po h 有关,用一阶泰勒展开式近似函数得到欧拉方法,其局部截断误差为一阶泰勒余项2()o h ,故是一阶方法,完全类似地若用p 阶泰勒展开式2'''()11()()()......()()2!!pp p n n n n n h h y y x hy x y x y x O h p ++=+++++ 进行离散化,所得计算公式必为p 阶方法,式中'''''()(,),()(,)(,)(,)....x y x f x y y x f x y f x y f x y ==++由此,我们能够想到,通过提高泰勒展开式的阶数,可以得到高精度的数值方法,从理论上讲,只要微分方程的解()y x 充分光滑,泰勒展开方法可以构造任意的有限阶的计算公式,但事实上,具体构造这种公式往往相当困难,因为符合函数(,())f x y x 的高阶导数常常是很烦琐的,因此,泰勒展开方法一般不直接使用,但是我们可以间接使用泰勒展开方法,求得高精度的计算方法。
首先,我们对欧拉公式和改进欧拉公式的形式作进一步的分析。
如果将欧拉公式和改进的欧拉公式改写成如下的形式:欧拉公式{111(,)n n n n y y hK K f x y +==+改进的欧拉公式11211()22n n y y h K K +=++, 1(,)n n K f x y =,21(,)n n K f x h y hK =++。
这两组公式都是用函数(,)f x y 在某些点上的值的线性组合来计算1()n y x +的近似值1n y +,欧拉公式每前进一步,就计算一次(,)f x y 的值。
另一方面它是1()n y x +在n x 处的一阶泰勒展开式,因而是一阶方法。
改进的欧拉公式每前进一步,需要计算两次(,)f x y 的值。
四阶龙格库塔法微分方程求解发散下载提示:该文档是本店铺精心编制而成的,希望大家下载后,能够帮助大家解决实际问题。
文档下载后可定制修改,请根据实际需要进行调整和使用,谢谢!本店铺为大家提供各种类型的实用资料,如教育随笔、日记赏析、句子摘抄、古诗大全、经典美文、话题作文、工作总结、词语解析、文案摘录、其他资料等等,想了解不同资料格式和写法,敬请关注!Download tips: This document is carefully compiled by this editor. I hope that after you download it, it can help you solve practical problems. The document can be customized and modified after downloading, please adjust and use it according to actual needs, thank you! In addition, this shop provides you with various types of practical materials, such as educational essays, diary appreciation, sentence excerpts, ancient poems, classic articles, topic composition, work summary, word parsing, copy excerpts, other materials and so on, want to know different data formats and writing methods, please pay attention!四阶龙格库塔法微分方程求解发散引言微分方程是数学中一种重要的工具,用于描述自然界中许多现象的变化规律。
四阶龙格库塔法求解微分方程组四阶龙格-库塔法(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))然后,我们通过逐步迭代的方式来求解方程组。
四阶龙格库塔法原理四阶龙格库塔法是一种常用的数值解微分方程的方法,它在工程、物理、计算机科学等领域都有着广泛的应用。
四阶龙格库塔法的原理相对简单,但能够较为准确地求解微分方程,因此受到了广泛的关注和应用。
首先,我们来了解一下什么是微分方程。
微分方程是描述自然界中变化规律的数学模型,它包含了未知函数及其导数的方程。
在实际问题中,我们常常需要求解微分方程来得到系统的演化规律。
然而,有些微分方程并没有解析解,这时就需要借助数值方法来求解。
四阶龙格库塔法就是其中一种常用的数值方法。
四阶龙格库塔法的原理可以简单概括为以下几个步骤:首先,我们需要将微分方程化为一阶方程组的形式。
这是因为四阶龙格库塔法是针对一阶方程组的数值解法。
其次,我们需要选择一个合适的步长h。
步长的选择对于数值解的精度有着重要的影响,通常需要通过实验来确定一个合适的步长。
然后,我们可以通过迭代的方式来逐步求解微分方程。
四阶龙格库塔法通过多次迭代,利用当前点的斜率来估计下一个点的值,从而逐步逼近真实解。
最后,我们可以通过迭代得到的数值解来近似地描述微分方程的解。
四阶龙格库塔法在一定条件下能够较为准确地求解微分方程,因此在实际应用中有着广泛的价值。
四阶龙格库塔法的原理相对简单,但是需要注意的是,在具体应用时需要考虑步长的选择、迭代次数的确定等参数,以及对于不同类型的微分方程可能需要进行适当的调整。
此外,四阶龙格库塔法也有其局限性,对于某些特殊类型的微分方程可能并不适用,因此在具体应用时需要综合考虑。
总的来说,四阶龙格库塔法是一种常用的数值解微分方程的方法,它通过迭代的方式逐步逼近真实解,具有一定的精度和稳定性。
在实际应用中,我们需要根据具体问题的特点来选择合适的数值方法,并注意参数的选择和调整,以获得准确的数值解。
四阶龙格库塔法在工程、物理、计算机科学等领域都有着广泛的应用,对于研究和解决实际问题具有重要的意义。
数值计算:四阶龙格-库塔法for⼆阶微分⽅程引⾔考虑存在以下⼆阶偏微分⽅程f2⋅¨X(t)+f1⋅˙X(t)+f0⋅X(t)=F(t)如何使⽤四阶龙格-库塔法求解该微分⽅程?⼀阶微分⽅程的解法⾸先回顾下对于⼀阶微分⽅程的解法,现在有以下⼀阶常系数⾮齐次微分⽅程f1⋅˙X(t)+f0⋅X(t)=F(t)⽅程可以写作˙X(t)=F(t)−f0⋅X(t)f1X n+1−X ndt=F(t)−f0⋅X(t)f1X n+1=X n+dt⋅F(t)−f0⋅X nf1=X n+dt⋅f(t,X n)其中dt为时间步长。
四阶龙格-库塔法公式如下:X n+1=X n+dt6k1+2k2+2k3+k4k1=f t,X nk2=f t+dt2,Xn+dt2⋅k1k3=f t+dt2,Xn+dt2⋅k2k4=f t+dt,X n+dt⋅k3利⽤以上公式即可显式推进求解⼆阶微分⽅程的解法现在考虑⼆阶常系数⾮齐次微分⽅程f2⋅¨X(t)+f1⋅˙X(t)+f0⋅X(t)=F(t)⽅程中出现⼆次导数项¨X(t),难以处理,所以考虑将⼆次导数项¨X(t)转化为⼀次导数项,⾄此,我们引⼊a(t)=X(t)b(t)=˙a(t)=˙X(t)原⽅程可以写成⼀阶偏微分⽅程组{()()()()(){f 2⋅˙b (t )+f 1⋅b (t )+f 0⋅a (t )=F (t )b (t )=˙a(t )仿照以上⼀阶微分⽅程的RK4解法˙a(t )=b (t )˙b (t )=F (t )−(f 1⋅b (t )+f 0⋅a (t ))f 2˙a n +1=a n +dt ⋅b (t )=a n +dt ⋅f (t ,a n ,b n )˙bn +1=b n +dt ⋅F (t )−(f 1⋅b (t )+f 0⋅a (t ))f 2=b n +dt ⋅g (t ,a n ,b n )˙a n +1=a n +dt ⋅f (t ,a n ,b n )˙bn +1=b n +dt ⋅g (t ,a n ,b n )化简⾄此,便可以仿照之前的RK4⽅法进⾏求解,具体公式为a n +1b n +1=a n +dt6⋅(k 1a +2k 2a +2k 3a +k 4a )b n +dt6⋅(k 1b +2k 2b +2k 3b +k 4b )其中,各个系数求解公式为:k 1a k 1b=f (t ,a n ,b n )g (t ,a n ,b n )k 2a k 2b =f (t +dt2,a n +dt2⋅k 1a ,b n +dt2⋅k 1b )g (t +dt2,a n +dt2⋅k 1a ,b n +dt2⋅k 1b )k 3a k 3b=f (t +dt2,a n +dt2⋅k 2a ,b n +dt2⋅k 2b )g (t +dt2,a n +dt2⋅k 2a ,b n +dt2⋅k 2b )k 4a k 4b=f (t +dt ,a n +dt ⋅k 3a ,b n +dt ⋅k 3b )g (t +dt ,a n +dt ⋅k 3a ,b n +dt ⋅k 3b )⾄此已经完成使⽤四阶龙格-库塔法⼆阶微分⽅程求解过程的梳理实例求解,以下偏微分⽅程:\begin{align} 1 \cdot \ddot{X(t)}+0 \cdot \dot{X(t)} +0 \cdot {X(t)} =cos(t)\\ \dot{X(0)}=0\\ {X(0)} =0 \end{align}理论解:\begin{align} {X(t)} =-cos(t)+1 \end{align}#include"RK_ODE.h"#include<iostream>using namespace std;#define M 1{{{{[][]{[][][][][][][][]MatrixXd F2(double t) {MatrixXd res = MatrixXd::Identity(M, M);return res;}MatrixXd F1(double t) {MatrixXd res = MatrixXd::Zero(M, M);return res;}MatrixXd F0(double t) {MatrixXd res = MatrixXd::Zero(M, M);return res;}MatrixXd F(double t) {MatrixXd res = MatrixXd::Ones(M, 1);res(0, 0) = cos(t);return res;}int main() {RK_ODE ode(M, 0.1);ode.X.fill(0.);ode.dX.fill(0.);for (int i = 0; i < 1000; i++) {ode.Solve(F2, F1, F0, F);ode.WriteMotionAndForce("results1.txt", 0); }//ode.SaveSnapshoot("snapshoot.json");cout << ode._Mode << endl;}Processing math: 88%。
一、介绍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中的四阶龙格库塔方法来求解微分方程组。
m a t l a b四阶龙格-库塔法求微分方程Matlab 实现四阶龙格-库塔发求解微分方程从理论上讲,只要函数在某区间上充分光滑,那么它可以展开为泰勒级数,因此在该区间上的函数值可用各阶导数值近似地表示出来,反之其各阶导数值也可用某些函数值的线性组合近似地表示出来。
龙格-库塔法就是将待求函数)(t y 展开为泰勒级数,并用方程函数),(y f t 近似其各阶导数,从而迭代得到)(t y 的数值解。
具体来说,四阶龙格-库塔迭代公式为)22(6143211k k k k h n n ++++=+y y ),(1n n t k y f =)2/,2/(12hk h t k n n ++=y f)2/,2/(23hk h t k n n ++=y f),(33hk h t k n n ++=y f实验内容:已知二阶系统21x x= ,u x x x 5.02.04.0212+--= ,0)0()0(21==x x ,u 为单位阶跃信号。
用四阶龙格-库塔法求数值解。
分析步长对结果的影响。
实验总结:实验报告要求简要的说明实验原理;简明扼要地总结实验内容;编制m 文件,并给出运行结果。
报告格式请按实验报告模板编写。
进入matlab ,Step1:choose way1 or way2way1):可以选择直接加载M 文件(函数M 文件)。
way2):点击new ——function ,先将shier (函数1文本文件)复制运行;点击new——function,再将RK(函数2文本文件)运行;点击new——function,再将finiRK(函数3文本文件)运行;Step2:回到command页面输入下面四句。
[t,k]=finiRK45([0;0],150);%迭代150次,步长=20/150[t1 k1]=ode45(@shier,[0 -10],[0 0]);%调用matlab自带四阶龙格-库塔,对比结果[t2 k2]=ode45(@shier,[0 10],[0 0]);plot(t,k(1,:),'-',t1,k1(:,1),'*',t2,k2(:,1),'^')%在图形上表示出来补充:改变步长影响数据的准确性。
经典四阶龙格库塔法解一阶微分方程组1.经典四阶龙格库塔法解一阶微分方程组1.1运用四阶龙格库塔法解一阶微分方程组算法分析(1-1),(1-2)(1-3)(1-4)(1-5)(1-6)(1-7)(1-8)(1-9)(1-10)经过循环计算由推得……每个龙格-库塔方法都是由一个合适的泰勒方法推导而来,使得其最终全局误差为,一种折中方法是每次进行若干次函数求值,从而省去高阶导数计算。
4阶龙格-库塔方法(RK4)是最常用的,它适用于一般的应用,因为它非常精准,稳定,且易于编程。
1.2经典四阶龙格库塔法解一阶微分方程流程图图1-1 经典四阶龙格库塔法解一阶微分方程流程图1.3经典四阶龙格库塔法解一阶微分方程程序代码:#include#includeusing namespace std;void RK4( double (*f)(double t,double x, double y),double (*g)(double t,double x, double y) ,double initial[3], double resu[3],double h){double f1,f2,f3,f4,g1,g2,g3,g4,t0,x0,y0,x1,y1;t0=initial[0];x0=initial[1];y0=initial[2];f1=f(t0,x0,y0); g1=g(t0,x0,y0);f2=f(t0+h/2, x0+h*f1/2,y0+h*g1/2); g2=g(t0+h/2, x0+h*f1/2,y0+h*g1/2);f3=f(t0+h/2, x0+h*f2/2,y0+h*g2/2); g3=g(t0+h/2, x0+h*f2/2,y0+h*g2/2);f4=f(t0+h, x0+h*f3,y0+h*g3); g4=g(t0+h, x0+h*f3,y0+h*g3); x1=x0+h*(f1+2*f2+2*f3+f4)/6; y1=y0+h*(g1+2*g2+2*g3+g4)/6; resu[0]=t0+h;resu[1]=x1;resu[2]=y1;}int main(){double f(double t,double x, double y);double g(double t,double x, double y);double initial[3],resu[3];double a,b,H;double t,step;int i;cout<<"输入所求微分方程组的初值t0,x0,y0:";cin>>initial[0]>>initial[1]>>initial[2];cout<<"输入所求微分方程组的微分区间[a,b]:";cin>>a>>b;cout<<"输入所求微分方程组所分解子区间的个数step:";cin>>step;cout<<setiosflags(ios::right)<<setiosflags(ios::fixed)<</seti osflags(ios::right)<<setiosflags(ios::fixed)<cout<<initial[0]<<setw(18)<<initial[1]<<setw(18)<<initial[2]<<endl;< bdsfid="133"p=""></setw(18)<<initial[1]<<setw(18)<<initial[2]<<endl;<> for(i=0;i<step;i++)< bdsfid="135" p=""></step;i++)<>{ RK4( f,g ,initial, resu,H);cout<<resu[0]<<setw(20)<<resu[1]<<setw(20)<<resu[2]< <endl;< bdsfid="138" p=""></resu[0]<<setw(20)<<resu[1]<<setw(20)<<resu[2]<<en dl;<>initial[0]=resu[0];initial[1]=resu[1];initial[2]=resu[2];}return(0);}double f(double t,double x, double y){double dx;dx=x+2*y;return(dx);}double g(double t,double x, double y)double dy;dy=3*x+2*y;return(dy);}1.4经典四阶龙格库塔法解一阶微分方程程序调试结果图示:应用所编写程序计算所给例题:其中初值为求解区间为[0,0.2]。
四阶龙格库塔法解微
分方程
四阶龙格库塔法解微分方程
一、四阶龙格库塔法解一阶微分方程
①一阶微分方程:
cos y
t & ,初始值y(0)=0,求解区间[0 10]。
MATLAB 程序:
%%%%%%%%%%% 四阶龙哥库塔法解一阶微分方程
%%%%%%%%%%% y'=cost
%%%%%%%%%%% y(0)=0, 0≤t ≤10,h=0.01
%%%%%%%%%%% y=sint
h=0.01;
hf=10;
t=0:h:hf;
y=zeros(1,length(t));
y(1)=0;
F=@(t,y)(cos(t));
for i=1:(length(t)-1)
k1=F(t(i),y(i));
k2=F(t(i)+h/2,y(i)+k1*h/2);
k3=F(t(i)+h/2,y(i)+k2*h/2);
k4=F(t(i)+h,y(i)+k3*h);
y(i+1)=y(i)+1/6*(k1+2*k2+2*k3+k4)*h;
end
subplot(211)
plot(t,y,'-')
xlabel('t');
ylabel('y');
title('Approximation');
span=[0,10];
p=y(1);
[t1,y1]=ode45(F,span,p);
subplot(212)
plot(t1,y1)
xlabel('t');
ylabel('y');
title('Exact');
图1
②一阶微分方程:()22*/x
t x x t =-&& ,初始值x(1)=2,求解区间[1
3]。
MATLAB 程序: %%%%%%%%%%% 四阶龙哥库塔法解微分方程
%%%%%%%%%%% x'(t)=(t*x-x^2)/t^2
%%%%%%%%%%% x(1)=2, 1≤t ≤3, h=1/128
%%%%%%%%%%% 精确解:x(t)=t/(0.5+lnt)
h=1/128; %%%%% 步长
tf=3;
t=1:h:tf;
x=zeros(1,length(t));
x(1)=2; %%%%% 初始值
F_tx=@(t,x)(t.*x-x.^2)./t.^2;
for i=1:(length(t)-1)
k_1=F_tx(t(i),x(i));
k_2=F_tx(t(i)+0.5*h,x(i)+0.5*h*k_1);
k_3=F_tx((t(i)+0.5*h),(x(i)+0.5*h*k_2)); k_4=F_tx((t(i)+h),(x(i)+k_3*h));
x(i+1)=x(i)+(1/6)*(k_1+2*k_2+2*k_3+k_4)*h; end
subplot(211)
plot(t,x,'-');
xlabel('t');
ylabel('x');
legend('Approximation');
%%%%%%%%%%%%%%%%%%%%%%%%%%%% ode45求精确解
t0=t(1);x0=x(1);
xspan=[t0 tf];
[x_ode45,y_ode45]=ode45(F_tx,xspan,x0);
subplot(212)
plot(x_ode45,y_ode45,'--');
xlabel('t');
ylabel('x');
legend('Exact');
图2
二、四阶龙格库塔法解二阶微分方程
①二阶微分方程:cos y t &
& ,初始值y(0)=0,y'(0)=-1,求解区间[0 10]。
MATLAB 程序:
%%%%%%%%% 龙格库塔欧拉方法解二阶微分方程
%%%%%%%%% y''=cost
%%%%%%%%% y'(0)=-1, y(0)=0
%%%%%%%%% 转换为一阶微分方程
h=0.01;
ht=10;
t=0:h:ht;
%%%%%%%%% y'=z1, z1'=y''=cost
y=zeros(1,length(t));
z=zeros(1,length(t));
y(1)=0;
z(1)=-1;
f=@(t,y,z)(z);
g=@(t,y,z)(sin(t));
for i=1:(length(t)-1)
K1=f(t(i),y(i),z(i));
L1=g(t(i),y(i),z(i));
K2=f(t(i)+h/2,y(i)+h/2*K1,z(i)+h/2*L1);
L2=g(t(i)+h/2,y(i)+h/2*K1,z(i)+h/2*L1);
K3=f(t(i)+h/2,y(i)+h/2*K2,z(i)+h/2*L2);
L3=g(t(i)+h/2,y(i)+h/2*K2,z(i)+h/2*L2);
K4=f(t(i)+h,y(i)+h*K3,z(i)+h*L3);
L4=g(t(i)+h,y(i)+h*K3,z(i)+h*L3);
y(i+1)=y(i)+h/6*(K1+2*K2+2*K3+K4);
z(i+1)=z(i)+h/6*(L1+2*L2+2*L3+L4);
end
plot(t,y)
xlabel('t');
ylabel('y');
title('y''=sin(t)');
legend('y''=sint');
图3
②二阶微分方程:7.35499*cos y x &
& ,初始值y(1)=0,y'(1)=0,求解区间[0 25]。
MATLAB 程序:
%%%%%%%%% 龙格库塔欧拉方法解二阶微分方程
%%%%%%%%% y''=7.35499*cosx
%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) %周期
function w=rightf_sys1(x,y,z)
w=7.35499*cos(y);
end
图4
注:四阶龙格库塔法求解二阶微分方程时,只需把二阶微分方程转化为一阶微分方程,再进行相关求解即可。