计算方法上机作业——龙格库塔法matlab程序
四阶龙格-库塔法求解常微分方程的初值问题-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。
内弹道 龙格库塔 计算 matlab
内弹道是指射程较短的导弹或火箭弹在飞行过程中受到大气阻力和重力等作用的飞行轨迹。
内弹道理论研究的是导弹或火箭弹在发射后到离开大气层再进入大气层末时的飞行过程。
内弹道包括导弹或火箭弹在发射后的加速、稳定、制导、飞行以及飞行过程中的动力学性能仿真等诸多内容。
内弹道有着复杂的飞行特性和动力学方程，在实际工程中需要进行准确的计算和仿真。
内弹道的计算中，龙格库塔（Runge-Kutta）法是一种常用的数值积分方法，在求解微分方程等领域有着广泛的应用。
龙格库塔法是由数学家奥特翁格（C. W. Runge）和马丁庫塔（M. W. J. Kutta）于1900年提出的，用于求解常微分方程初值问题，其优点是精度较高，适用范围广。
在内弹道计算中，可以利用龙格库塔法对导弹或火箭弹的飞行轨迹进行数值模拟和计算，得到较为准确的飞行轨迹数据。
在实际工程中，为了方便进行内弹道的计算，可以使用Matlab等数学建模和仿真软件。
Matlab是一种常用的科学计算软件，具有强大的数值计算和仿真功能，可以用于内弹道计算中的龙格库塔法数值模拟。
在Matlab中，可以编写相应的程序，利用龙格库塔法对导弹或火箭弹的飞行过程进行仿真和计算，得到准确的飞行轨迹和动力学性能数据。
内弹道计算是导弹或火箭弹研究设计中的重要内容，龙格库塔法是一种常用的数值积分方法，Matlab是一种常用的科学计算软件，它们的应用能够有效地进行内弹道的计算和仿真，为导弹或火箭弹的研制提供重要的技术支持。
随着技术的不断发展，内弹道计算已经成为导弹或火箭弹研究设计中不可或缺的一部分。
在内弹道计算中，龙格库塔法是一种常用的数值积分方法，可以对导弹或火箭弹的飞行轨迹进行数值模拟和计算，提供准确的飞行轨迹数据。
而Matlab作为一种强大的科学计算软件，对于内弹道的计算和仿真也有着重要的应用价值。
在实际工程中，使用Matlab编写程序，利用龙格库塔法对导弹或火箭弹的飞行轨迹进行数值模拟和计算，将为导弹或火箭弹的研制提供重要的技术支持。
matlab龙格库塔法程序,给出实例
一、介绍龙格库塔法龙格库塔法（Runge-Kutta method）是一种数值计算方法，用于求解常微分方程的数值解。
它通过多步迭代的方式逼近微分方程的解，并且具有较高的精度和稳定性。
二、龙格库塔法的原理龙格库塔法采用迭代的方式来逼近微分方程的解。
在每一步迭代中，计算出当前时刻的斜率，然后根据这个斜率来求解下一个时刻的值。
通过多步迭代，可以得到微分方程的数值解。
三、龙格库塔法的公式龙格库塔法可以表示为以下形式：k1 = f(tn, yn)k2 = f(tn + h/2, yn + h/2 * k1)k3 = f(tn + h/2, yn + h/2 * k2)k4 = f(tn + h, yn + h * k3)yn+1 = yn + h/6 * (k1 + 2k2 + 2k3 + k4)其中，k1、k2、k3、k4为斜率，h为步长，tn为当前时刻，yn为当前时刻的解，yn+1为下一个时刻的解。
四、使用matlab实现龙格库塔法在MATLAB中，可以通过编写函数来实现龙格库塔法。
下面是一个用MATLAB实现龙格库塔法的简单例子：```matlabfunction [t, y] = runge_kutta(f, tspan, y0, h)t0 = tspan(1);tf = tspan(2);t = t0:h:tf;n = length(t);y = zeros(1, n);y(1) = y0;for i = 1:n-1k1 = f(t(i), y(i));k2 = f(t(i) + h/2, y(i) + h/2 * k1);k3 = f(t(i) + h/2, y(i) + h/2 * k2);k4 = f(t(i) + h, y(i) + h * k3);y(i+1) = y(i) + h/6 * (k1 + 2*k2 + 2*k3 + k4);endend```以上就是一个简单的MATLAB函数，可以利用该函数求解给定的微分方程。
微分方程的数值解法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-欧拉方法和龙格库塔方法的小实例
题一：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 应用使用Euler 和Rungkutta 方法解臂状摆的能量方程背景单摆是常见的物理模型，为了得到摆角θ的关于时间的函数，来描述单摆运动。
由角动量定理我们知道εJ M =化简得到 0sin 22=+θθlg dt d 在一般的应用和计算中，只考虑摆角在5度以内的小摆动，因为可以吧简化为θ，这样比较容易解。
实际上这是一个解二阶常微分方程的问题。
在这里的单摆是一种特别的单摆，具有均匀的质量M 分布在长为2的臂状摆上，使用能量法建立方程 WT = h mg ?=2J 21ω 化简得到θθcos35499.722=dtd 重力加速度取9.806651使用欧拉法令dxdy z =，这样降阶就把二阶常微分方程转化为一阶微分方程组，再利用向前Euler 方法数值求解。
y(i+1)=y(i)+h*z(i);z(i+1)=z(i)+h*7.35499*cos(y(i));y(0)=0z(0)=0精度随着h 的减小而更高，因为向前欧拉方法的整体截断误差与h 同阶，（因为是用了泰勒公式）所以欧拉方法的稳定区域并不大。
2.RK4-四阶龙格库塔方法使用四级四阶经典显式Rungkutta 公式稳定性很好，RK4法是四阶方法，每步的误差是h5阶，而总积累误差为h4阶。
所以比欧拉稳定。
运行第三个程序：在一幅图中显示欧拉法和RK4法，随着截断误差的积累，欧拉法产生了较大的误差h=0.01h=0.0001，仍然是开始较为稳定，逐渐误差变大总结：RK4是很好的方法，很稳定，而且四阶是很常用的方法，因为到五阶的时候精度并没有相应提升。
通过这两种方法计算出角度峰值y=3.141593，周期是1.777510。
三个程序欧拉法clear;clch=0.00001;a=0;b=25;x=a:h:b;y(1)=0;z(1)=0;for i=1:length(x)-1 % 欧拉y(i+1)=y(i)+h*z(i);z(i+1)=z(i)+h*7.35499*cos(y(i));endplot(x,y,'r*');xlabel('时间');ylabel('角度');A=[x,y];%y(find(y==max(y)))%Num=(find(y==max(y)))[y,T]=max(y);fprintf('角度峰值等于%d',y) %角度的峰值也就是πfprintf('\n')fprintf('周期等于%d',T*h) %周期legend('欧拉');龙格库塔方法先定义函数rightf_sys1.mfunction w=rightf_sys1(x,y,z)w=7.35499*cos(y);clear;clc;%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)-1K1=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 L4RK_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);endplot(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) %周期两个方法在一起对比使用跟上一个相同的函数rightf_sys1.mclear;clc; %清屏h=0.0001;a=0;b=25;x=a:h:b;Euler_y(1)=0;Euler_z(1)=0; %欧拉的初值RK_y(1)=0;RK_z(1)=0; %龙格库塔初值for i=1:length(x)-1%先是欧拉法Euler_y(i+1)=Euler_y(i)+h*Euler_z(i);Euler_z(i+1)=Euler_z(i)+h*7.35499*cos(Euler_y(i));%龙格库塔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);% K2 and L2K3=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 );% K3 and L3K4=RK_z(i)+h*L3;L4=rightf_sys1(x(i)+h,RK_y(i)+h*K3,RK_z(i)+h*L3); % K4 and L4RK_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);endplot(x,Euler_y,'r-',x,RK_y,'b-');[y,T]=max(RK_y);fprintf('角度峰值等于%d',y) %角度的峰值也就是πfpri ntf('\n') fprintf('周期等于%d',T*h) %周期xlabel('时间');ylabel('角度');legend('欧拉','RK4');。
常微分方程组的四阶RUNGEKUTTA龙格库塔法MATLAB实现
常微分方程组的四阶RUNGEKUTTA龙格库塔法MATLAB实现欢迎使用 Python 版的实现！数值解常微分方程组的四阶 Runge-Kutta 方法（也被称为 RK4 方法）可以通过迭代的方式逐步逼近精确解。
对于一阶常微分方程组：dy/dt = f(t, y)我们可以通过下面的公式来计算下一个时间步长上的近似解：y(n+1)=y(n)+(1/6)(k1+2k2+2k3+k4)其中k1=h*f(t(n),y(n))k2=h*f(t(n)+h/2,y(n)+k1/2)k3=h*f(t(n)+h/2,y(n)+k2/2)k4=h*f(t(n)+h,y(n)+k3)这里，h代表时间步长，t(n)代表当前时间，y(n)代表当前解。
f(t,y)是给定的方程组。
对于四阶 Runge-Kutta 方法的 MATLAB 实现，可以按照以下步骤进行。
1.首先，定义需要求解的常微分方程组。
function dydt = equations(t, y)dydt = zeros(2, 1);dydt(1) = y(2); % 根据方程组的具体形式修改dydt(2) = -y(1); % 根据方程组的具体形式修改end2.定义RK4方法的求解函数。
function [t, y] = rk4_solver(equations, tspan, y0, h) t = tspan(1):h:tspan(2);y = zeros(length(y0), length(t));y(:,1)=y0;for i = 1:length(t)-1k1 = h * equations(t(i), y(:, i));k2 = h * equations(t(i) + h/2, y(:, i) + k1/2);k3 = h * equations(t(i) + h/2, y(:, i) + k2/2);k4 = h * equations(t(i) + h, y(:, i) + k3);y(:,i+1)=y(:,i)+(k1+2*k2+2*k3+k4)/6;endend3. 调用 rk4_solver 函数求解常微分方程组。
计算方法上机作业——龙格库塔法
1 yi 1 yi 6 K1 2 K 2 2 K 3 K 4 K hf x , y i i 1 K1 h K 2 hf xi , yi 2 2 K2 h K 3 hf xi 2 , yi 2 K hf x h, y K 4 i i 3 1 y j ,i 1 y j ,i 6 K j1 2 K j 2 2 K j 3 K j 4 K hf x ; y , y , , y i 1i 2i ni j1 K K11 K h , y2i 21 , , yni n1 K j 2 hf xi ; y1i 2 2 2 2 Kn2 K12 K 22 h K j 3 hf xi 2 ; y1i 2 , y2i 2 , , yni 2 K hf x h; y K , y K , , y K j4 i 1i 13 2i 23 ni n3
四阶龙格库塔法求解常微分方程的初值问题18图19标准四级四阶龙格库塔法程序框图42程序使用说明标准四级四阶龙格库塔法求解常微分方程初值问题的matlab程序见附录程序可以用来求解一阶常微分方程组和高阶常微分方程的初值问题运行程序按照命令窗口的提示输入相关变量直至得到结果
计算方法上机报告
4 四阶龙格－库塔法求解常微分方程的初值问题
4.1 算法原理及程序框图 一阶常微分方程（组）和高阶常微分方程的初值问题最终都可以转化为一阶常微 分方程组的初值问题，其向量形式为式(17)。
y x f x, y x , y a y0
a xb
(17)
式(17)在形式上与单个微分方程的初值问题完全相同， 只是将数量函数变成了向量 函数。 因此， 标准四级四阶龙格—库塔法的向量形式和分量形式分别为式(18)和式(19)。
matlab ppt_13_龙格-库塔法
这就是改进的欧拉公式，通常写成 h y = y + ( K1 + K2 ) i i+1 2 K1 = f ( x i , y i ) K2 = f (xi + h, yi + hK1) 由此可见，为了提高精度，可以多取几点的斜率值作加权平均当作平均斜 率kave，这就是龙格库塔法基本思想. 一般的显式龙格库塔方法可以写成
2. 二阶龙格――库塔法
在改进的欧拉公式中将xi+1替换为区间中的任意一点 xi+p = xi + ph (0 < p ≤ 1)
将它的斜率记为K2。利用欧拉公式得y (xi+p)的预报值 y i+p = yi + phK1 求得 K2 = f (xi+p, y i+p) 平均斜率则取 Keve ≈ λ1K1 + λ2K2 其中λ1, λ2是待定的系数。
(h)
y (xi+1) − yi+1 由此可得下列事后估计式
2 y (xi+1) − yi+1 ≈
(h)
≈
1 16
(h)
1 (h (h) 2) (yi+1 − yi+1) 15
这样，可以通过检查步长折半前后两次计算结果的偏差 ∆=
(h 2) |yi+1
− yi+1|
(j )
来判断所选取的步长是否合适。具体可分以下两种情况： 1. 对于给定的精度 ,如果∆ > ，就反复将步长折半计算，直至∆ < 为 (h 2) 止，这时取步长折半后的“新值”yi+1 作为结果。 2. 如果∆ < ，就反复将步长加倍，直至∆ > 为止，这时取步长加倍前 的“旧值”作为结果。 这种通过步长折半或加倍的计算的方法就叫变步长方法。表面上看，为了 选择步长，每一步的计算量似乎增加了，但从总体上考虑往往是合算的。
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`，该向量包含时间范围的起始和终止值。
matlab龙格库塔法解动力学方程
matlab龙格库塔法解动力学方程
MATLAB提供了几个函数来实现龙格库塔法解动力学方程。
以下是一种基本方法的示例：
1. 定义动力学方程的函数
```matlab
function dy = dynamics(t, y)
g = 9.8; % 重力加速度
m = 1; % 物体质量
k = 0.1; % 阻尼系数
dy = zeros(2,1);
dy(1) = y(2);
dy(2) = -g - k/m * y(2);
end
```
2. 定义初始条件和时间范围
```matlab
y0 = [0; 0]; % 初始位置和速度
tspan = [0 5]; % 时间范围
```
3. 调用 ode45 函数来求解动力学方程
```matlab
[t, y] = ode45(@dynamics, tspan, y0);
```
4. 绘制位置和速度随时间变化的图形
```matlab
figure;
subplot(2,1,1);
plot(t, y(:,1));
xlabel('时间');
ylabel('位置');
subplot(2,1,2);
plot(t, y(:,2));
xlabel('时间');
ylabel('速度');
```
这个例子使用了ode45函数来解决动力学方程。
可以使用不同的龙格库塔法，如ode23或ode113，具体取决于问题的性质。
龙格-库塔法MATLAB
1. matlab 新建.m文件，编写龙格-库塔法求解函数function [x,y]=runge_kutta1(ufunc,y0,h,a,b)%参数表顺序依次是微分方程组的函数名称，初始值向量，步长，时间起点，时间终点（参数形式参考了ode45函数）n=floor((b-a)/h); %求步数x(1)=a;%时间起点y(:,1)=y0;%赋初值，可以是向量，但是要注意维数for ii=1:nx(ii+1)=x(ii)+h;k1=ufunc(x(ii),y(:,ii));k2=ufunc(x(ii)+h/2,y(:,ii)+h*k1/2);k3=ufunc(x(ii)+h/2,y(:,ii)+h*k2/2);k4=ufunc(x(ii)+h,y(:,ii)+h*k3);y(:,ii+1)=y(:,ii)+h*(k1+2*k2+2*k3+k4)/6;%按照龙格库塔方法进行数值求解end2.另外再新建一个.,m文件，定义要求解的常微分方程函数function dx=fun1(t,x)dx =zeros(2,1);%初始化列向量dx(1) =0.08574*x(2)-1.8874*x(1)-20.17;dx(2) =1.8874*x(1)-0.08574*x(2)+115.16;3,再新建一个.m文件，利用龙格-库塔方法求解常微分方程[T1,F1]=runge_kutta1(@fun1,[46.30 1296 ],1,0,20); %求解步骤2定义的fun1常微分方程，@fun1是调用其函数指针，从0到20，间隔为1subplot(122)plot(T1,F1)%自编的龙格库塔函数效果title('自编的龙格库塔函数')grid运行步骤3文件即可得到结果，F1为估测值或者可以调用matlab自带函数ode45求解方法如下：[T,F] = ode45(@fun1,0:1:20,[17.64 37800 ]); subplot(121)plot(T,F)%Matlab自带的ode45函数效果title('ode45函数效果')grid。
龙格库塔方法及其matlab实现
资料范本本资料为word版本，可以直接编辑和打印，感谢您的下载龙格库塔方法及其matlab实现地点：__________________时间：__________________说明：本资料适用于约定双方经过谈判，协商而共同承认，共同遵守的责任与义务，仅供参考，文档可直接下载或修改，不需要的部分可直接删除，使用时请详细阅读内容龙格-库塔方法及其matlab实现摘要：本文的目的数值求解微分方程精确解，通过龙格-库塔法，加以利用matlab为工具达到求解目的。
龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法，用于数值求解微分方程。
MatLab软件是由美国Mathworks公司推出的用于数值计算和图形处理的科学计算系统环境。
MatLab 是英文MATrix LABoratory（矩阵实验室）的缩写。
在MratLab环境下，用户可以集成地进行程序设计、数值计算、图形绘制、输入输出、文件管理等各项操作。
关键词：龙格-库塔 matlab 微分方程前言1.1：知识背景龙格－库塔法（Runge-Kutta）是用于非线性常微分方程的解的重要的一类隐式或显式迭代法。
这些技术由数学家卡尔·龙格和马丁·威尔海姆·库塔于1900年左右发明。
通常所说的龙格库塔方法是相对四阶龙格库塔而言的，成为经典四阶龙格库塔法。
该方法具有精度高，收敛，稳定，计算过程中可以改变步长不需要计算高阶导数等优点，但是仍需计算在一些点上的值，比如四阶龙格-库塔法没计算一步需要计算四步，在实际运用中是有一定复杂性的。
Matlab是在20世纪七十年代后期的事：时任美国新墨西哥大学计算机科学系主任的Cleve Moler教授出于减轻学生编程负担的动机，为学生设计了一组调用LINPACK和EISPACK库程序的“通俗易用”的接口，此即用FORTRAN编写的萌芽状态的MATLAB。
经几年的校际流传，在Little的推动下，由Little、Moler、Steve Bangert合作，于1984年成立了MathWorks公司，并把MATLAB正式推向市场。
(完整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由图看出龙格库塔方法误差很小，具有很高的精度。
Matlab数值计算龙格库塔
clearclch=0.001;%前几行都是给变量赋初始值，是已知的x0=0;y0=0;t0=0;for i=1:0.1/h %循环，从起点1到终点0.1/h，每循环一次，i增大1，以下循环是为我军导弹服务的t=t0+h/2;%龙格公式里，每一个步长，分为起点、中点、终点，比如，t0是起点，t 是中点、t1是终点t1=t0+h;d1=450/(1+((90*t0*sin(0.3*pi)-y0)/(30+90*t0*cos(0.3*pi)-x0))^2)^(1/2);%龙格公式里坐标点（t0，x0）的导数,也是该点斜率d11=450/(1+((90*t0*sin(0.3*pi)-y0)\(30+90*t0*cos(0.3*pi)-x0))^2)^(1/2);%龙格公式里坐标点（t0，y0）的导数，也是该点斜率x=x0+h/2*d1;%龙格公式第一步，计算x的y=y0+h/2*d11;;%龙格公式第一步，计算yd2=450/(1+((90*t*sin(0.3*pi)-y)/(30+90*t*cos(0.3*pi)-x))^2)^(1/2);%龙格公式里坐标点（t，x）的导数d22=450/(1+((90*t*sin(0.3*pi)-y)\(30+90*t*cos(0.3*pi)-x))^2)^(1/2);%龙格公式里坐标点（t，y）的导数xx=x0+h/2*d2;%龙格公式第二步，计算x的yy=y0+h/2*d22;%龙格公式第二步，计算y的d3=450/(1+((90*t*sin(0.3*pi)-yy)/(30+90*t*cos(0.3*pi)-xx))^2)^(1/2);%龙格公式里坐标点（t，xx）的导数d33=450/((1+(90*t*sin(0.3*pi)-yy)\(30+90*t*cos(0.3*pi)-xx))^2)^(1/2);%龙格公式里坐标点（t，y）的导数x1x1=x0+h*d3;;%龙格公式第三步，计算x的y1y1=y0+h*d33;;%龙格公式第三步，计算y的d4=450/(1+((90*t1*sin(0.3*pi)-y1y1)/(30+90*t0*cos(0.3*pi)-x1x1))^2)^(1/2);%龙格公式里坐标点（t，x1x1）的导数d44=450/(1+((90*t1*sin(0.3*pi)-y1y1)\(30+90*t1*cos(0.3*pi)-x1x1))^2)^(1/2);%龙格公式里坐标点（t，y1y1）的导数x1(i)=x0+h/6*(d1+2*d2+2*d3+d4);%龙格公式第四步，计算x的y1(i)=y0+h/6*(d11+2*d22+2*d33+d44);%龙格公式第四步，计算y的t0=t1;%现在一个步长的龙格公式已结束，该到要进入下一个步长。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
计算方法上机报告 for mm=1:n fprintf('%9c%.f','y',mm) end for nn=1:N fprintf('\n%10.5f',x(nn)); for ll=1:n fprintf('%10.5f',y(ll,nn)); end if nn==N fprintf('\n'); end tlab 程序
附录 5 龙格—库塔法求解常微分方程初值问题的 matlab 程序
clear;clc; %======需要输入的=======% n = input('\n 请输入一阶微分方程组的个数或者高阶方程组的阶数 n：\n'); a = input('\n 请输入求解区间的下限 a:\n'); b = input('\n 请输入求解区间的上限 b:\n'); h = input('\n 请输入求解步长 h:\n'); fprintf('\n 请依次输入%.f 个方程的右端函数 fi(x,y(x)):\n',n) s = cell(n,1); for p=1:n fprintf('f%.f(x,y(x)) =',p) r = input(' ','s'); s{p} = r; end fprintf('\n 请依次输入%.f 个初值 yi(%.f):\n',n,a) y0 = zeros(n,1); for q=1:n fprintf('y%.f(%.f) =',q,a); y0(q) = input(' '); end %======需要输入的=======% N = (b-a)/h+1; x = a:h:b; K = zeros(n,4); y = zeros(n,N); y(:,1) = y0; for i=1:N-1 for a=1:n f = inline(char(s(a)),'x','y'); K(a,1) = h*f(x(i),y(:,i)); end for b=1:n f = inline(char(s(b)),'x','y'); K(b,2) = h*f(x(i)+h/2,y(:,i)+K(:,1)/2); end for c=1:n f = inline(char(s(c)),'x','y'); K(c,3) = h*f(x(i)+h/2,y(:,i)+K(:,2)/2); end for d=1:n f = inline(char(s(d)),'x','y'); K(d,4) = h*f(x(i)+h,y(:,i)+K(:,3)); end for j=1:n y(j,i+1) = y(j,i)+1/6*(K(j,1)+2*K(j,2)+2*K(j,3)+K(j,4)); end plot(x(i),y(1,i),'.k','markersize',15) hold on; end grid; fprintf('\n%10c','x')
29