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

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

一、四阶龙格库塔法解一阶微分方程

①一阶微分方程：

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

注：四阶龙格库塔法求解二阶微分方程时，只需把二阶微分方程转化为一阶微分方程，再进行相关求解即可。

龙格库塔解微分方程

《数值分析》课程实验报告 一、实验目的 1．掌握用MA TLAB求微分方程初值问题数值解的方法； 2．通过实例学习微分方程模型解决简化的实际问题； 3．了解龙格库塔方法的基本思想。 二、实验内容 用龙格库塔方法求下列微分方程初值问题的数值解 y’=x+y y(0)=1 0

end fprintf(‘结果矩阵，第一列为x(n),第二列~第五列为k1~k4，第六列为y（n+1）的结果') z %在命令框输入下列语句 %yy=inline('x+y'); %>> rk(yy,0,1,0.2,0,1) %将得到结果 四、实验小结 通过实验结果分析可知，龙格库塔方法求解常微分方程能获得比较好的数值解，龙格库塔方法的数值解的精度较高，方法比较简便易懂。

MATLAB改进欧拉法与四阶龙格-库塔求解一阶常微分方程

姓名：樊元君学号：02 日期： 一、实验目的 掌握MATLAB语言、C/C++语言编写计算程序的方法、掌握改进欧拉法与四阶龙格-库塔求解一阶常微分方程的初值问题。掌握使用MATLAB程序求解常微分方程问题的方法。 ： 二、实验内容 1、分别写出改进欧拉法与四阶龙格-库塔求解的算法，编写程序上机调试出结果，要求所编程序适用于任何一阶常微分方程的数值解问题，即能解决这一类问题，而不是某一个问题。 实验中以下列数据验证程序的正确性。 求，步长h=。 * 2、实验注意事项 的精确解为，通过调整步长，观察结果的精度的变化 ^ )

三、程序流程图： ●改进欧拉格式流程图： ~ |

●四阶龙格库塔流程图： ] 四、源程序： ●改进后欧拉格式程序源代码： function [] = GJOL(h,x0,y0,X,Y) format long h=input('h='); … x0=input('x0='); y0=input('y0='); disp('输入的范围是：'); X=input('X=');Y=input('Y='); n=round((Y-X)/h); \

i=1;x1=0;yp=0;yc=0; for i=1:1:n x1=x0+h; yp=y0+h*(-x0*(y0)^2);%yp=y0+h*(y0-2*x0/y0);% · yc=y0+h*(-x1*(yp)^2);%yc=y0+h*(yp-2*x1/yp);% y1=(yp+yc)/2; x0=x1;y0=y1; y=2/(1+x0^2);%y=sqrt(1+2*x0);% fprintf('结果=%.3f,%.8f,%.8f\n',x1,y1,y); : end end ●四阶龙格库塔程序源代码： function [] = LGKT(h,x0,y0,X,Y) 。 format long h=input('h='); x0=input('x0='); y0=input('y0='); disp('输入的范围是：'); " X=input('X=');Y=input('Y='); n=round((Y-X)/h); i=1;x1=0;k1=0;k2=0;k3=0;k4=0; for i=1:1:n ~ x1=x0+h; k1=-x0*y0^2;%k1=y0-2*x0/y0;% k2=(-(x0+h/2)*(y0+h/2*k1)^2);%k2=(y0+h/2*k1)-2*(x0+h/2)/(y0+h/2*k1);% k3=(-(x0+h/2)*(y0+h/2*k2)^2);%k3=(y0+h/2*k2)-2*(x0+h/2)/(y0+h/2*k2);% k4=(-(x1)*(y0+h*k3)^2);%k4=(y0+h*k3)-2*(x1)/(y0+h*k3);% … y1=y0+h/6*(k1+2*k2+2*k3+k4);%y1=y0+h/6*(k1+2*k2+2*k3+k4);% x0=x1;y0=y1; y=2/(1+x0^2);%y=sqrt(1+2*x0);% fprintf('结果=%.3f,%.7f,%.7f\n',x1,y1,y); end · end

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

matlab编的4阶龙格库塔法解微分方程的程序 2010-03-10 20:16 function varargout=saxplaxliu(varargin) clc,clear x0=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:n fprintf('%2d %4.4f %4.4f\n',i,x(i),y(i)); end function 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-1 K1=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); end y=y1; 结果： i x(i) y(i) 1 0.0000 1.0000 2 0.1000 0.9901 3 0.2000 0.9615 4 0.3000 0.9174 5 0.4000 0.8621 6 0.5000 0.8000 7 0.6000 0.7353 8 0.7000 0.6711 9 0.8000 0.6098 10 0.9000 0.5525 11 1.0000 0.5000 12 1.1000 0.4525 13 1.2000 0.4098

matlab 四阶龙格-库塔法求微分方程

Matlab 实现四阶龙格-库塔发求解微分方程 从理论上讲，只要函数在某区间上充分光滑，那么它可以展开为泰勒级数，因此在该区间上的函数值可用各阶导数值近似地表示出来，反之其各阶导数值也可用某些函数值的线性组合近似地表示出来。龙格-库塔法就是将待求函数)(t y 展开为泰勒级数，并用方程函数),(y f t 近似其各阶导数，从而迭代得到)(t y 的数值解。具体来说，四阶龙格-库塔迭代公式为 )22(6 143211k 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 way2 way1）： 可以选择直接加载M 文件（函数M 文件）。 way2）： 点击new ——function ，先将shier （函数1文本文件）复制运行; 点击new ——function ，再将RK （函数2文本文件）运行; 点击new ——function ，再将finiRK （函数3文本文件）运行;

经典四阶龙格库塔法解一阶微分方程组

1.经典四阶龙格库塔法解一阶微分方程组 1.1运用四阶龙格库塔法解一阶微分方程组算法分析 ),,(1k k k y x t f f =， )2,2,2(112g h y f h x h t f f k k k +++= )2 ,2,2(223g h y f h x h t f f k k k +++= ),,(334hg y hf x h t f f k k k +++= ),,(1k k k y x t g g = )2,2,2(112g h y f h x h t g g k k k +++= )2,2,2(223g h y f h x h t g g k k k +++= ),,(334hg y hf x h t g g k k k +++= ) 22(6 )22(6 43211 43211g g g g h y y f f f f h x x k k k k ++++=++++=++ 1k k t t h +=+ 经过循环计算由 推得 …… 每个龙格-库塔方法都是由一个合适的泰勒方法推导而来，使得其最终全局误差为() N O h ,一种折中方法是每次进行若干次函数求值，从而省去高阶导数计算。4阶龙格-库塔方法(RK4)是最常用的，它适用于一般的应用，因为它非常精 准，稳定，且易于编程。 000,,t x y ()()111222,,,,t x y t x y （1-1） （1-2） （1-3） （1-4） （1-5） （1-6） （1-7） （1-8） （1-9） （1-10）

1.2经典四阶龙格库塔法解一阶微分方程流程图 图1-1 经典四阶龙格库塔法解一阶微分方程流程图 1.3经典四阶龙格库塔法解一阶微分方程程序代码： #include #include using 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,

MATLAB龙格-库塔方法解微分方程

龙格-库塔方法是一种经典方法，具有很高的精度，它间接的利用了泰勒级数展开，避免了高阶偏导数的计算。此处以最为经典的四级四阶龙格-库塔方法为例，计算格式如下 ()()()11234121324 3226,,22,+22,n n n n n n n n n n h y y K K K K K f x y h h K f x y K h h K f x y K K f x h y hK +?=++++???=????=++? ???????=+? ?????=++? 1龙格-库塔法解一阶ODE 对于形如()()0, dy f x y a x b dx y a y ?=<≤???=? 的一阶ODE 初值问题，可以直接套用公式，如今可以借助计算机方便的进行计算，下面给出一个实例 ()2 0101dy x y x dx y y ?=-<≤???=? 取步长h=0.1 ，此处由数学知识可得该方程的精确解为y =。在这里利用MATLAB 编程，计算数值解并与精确解相比，代码如下： (1)写出微分方程，便于调用和修改 function val = odefun( x,y ) val = y-2*x/y; end (2)编写runge-kutta 方法的函数代码

function y = runge_kutta( h,x0,y0 ) k1 = odefun(x0,y0); k2 = odefun(x0+h/2,y0+h/2*k1); k3 = odefun(x0+h/2,y0+h/2*k2); k4 = odefun(x0+h,y0+h*k3); y = y0+h*(k1+2*k2+2*k3+k4)/6; end (3)编写主函数解微分方程,并观察数值解与精确解的差异clear all h = 0.1; x0 = 0; y0 = 1; x = 0.1:h:1; y(1) = runge_kutta(h,x0,y0); for k=1:length(x) x(k) = x0+k*h; y(k+1) = runge_kutta(h,x(k),y(k)); end z = sqrt(1+2*x); plot(x,y,’*’);

龙格库塔算法解微分方程组-c语言

This program is to solve the initial value problem of following system of differential equations: dx/dt=x+2*y,x(0)=0, dy/dt=2*x+y,y(0)=2, x and y are to be calculated ****************************************************************************/ #include<> #include<> #define steplength //步长?è可¨|根¨′据Y需¨¨要°a调ì??整; #define FuncNumber 2 //FuncNumber为a微?é分¤方¤程¨?的ì数oy目; void main() { float x[200],Yn[20][200],reachpoint;int i; x[0]=0;Yn[0][0]=0;Yn[1][0]=2; //初值|ì条?件t; reachpoint=; //所¨′求¨?点ì可¨|根¨′据Y需¨¨要°a调ì??整; void rightfunctions(float x ,float *Auxiliary,float *Func); void Runge_Kutta(float *x,float reachpoint, float(*Yn)[200]); Runge_Kutta(x ,reachpoint, Yn); printf("x "); for(i=0;i<=(reachpoint-x[0])/steplength;i++) printf("%f ",x[i]); printf("\nY1 "); for(i=0;i<=(reachpoint-x[0])/steplength;i++) printf("%f ",Yn[0][i]); printf("\nY2 "); for(i=0;i<=(reachpoint-x[0])/steplength;i++) printf("%f ",Yn[1][i]); getchar(); } void rightfunctions(float x ,float *Auxiliary,float *Func)//当ì?à右?¨°方¤程¨?改变à时o?à，ê需¨¨要°a改变à; { Func[0]=Auxiliary[0]+2*Auxiliary[1]; Func[1]=2*Auxiliary[0]+Auxiliary[1]; } void Runge_Kutta(float *x,float reachpoint, float(*Yn)[200]) { int i,j; float Func[FuncNumber],K[FuncNumber][4],Auxiliary[FuncNumber]; for(i=0;i<=(reachpoint-x[0])/steplength;i++) { for(j=0;j

龙格库塔算法解微分方程组 c语言

/*************************************************************************** This program is to solve the initial value problem of following system of differential equati ons: dx/dt=x+2*y,x(0)=0, dy/dt=2*x+y,y(0)=2, x(0.2) and y(0.2) are to be calculated ****************************************************************************/ #include #include #define steplength 0.1 //步?长?è可¨|根¨′据Y需¨¨要°a调ì??整?; #define FuncNumber 2 //FuncNumber为a微?é分¤?方¤?程¨?的ì?数oy目?; void main() { float x[200],Yn[20][200],reachpoint;i nt i; x[0]=0;Yn[0][0]=0;Yn[1][0]=2; //初?值|ì条??件t; reachpoint=0.2; //所¨′求¨?点ì?可¨|根¨′据Y需¨¨要°a调ì??整?; void rightfunctions(float x ,float *Auxiliary,float *Func); void R unge_Kutta(float *x,float reachpoint, float(*Yn)[200]); Runge_Kutta(x ,reachpoi nt, Yn); printf("x "); for(i=0;i<=(reachpoint-x[0])/steplength;i++) printf("%f ",x[i]); printf("\nY1 "); for(i=0;i<=(reachpoint-x[0])/steplength;i++) printf("%f ",Yn[0][i]); printf("\nY2 "); for(i=0;i<=(reachpoint-x[0])/steplength;i++) printf("%f ",Yn[1][i]); getchar(); } void rightfunctions(float x ,float *Auxiliary,float *Func)//当ì?à右?¨°方¤?程¨?改?变à?时o?à，ê?需¨¨要°a改?变à?; { Func[0]=Auxiliary[0]+2*Auxiliary[1]; Func[1]=2*Auxiliary[0]+Auxiliary[1]; } void R unge_Kutta(float *x,float reachpoint, float(*Yn)[200]) { int i,j; float Func[FuncNumber],K[FuncNumber][4],Auxiliary[FuncNumber]; for(i=0;i<=(reachpoint-x[0])/steplength;i++) { for(j=0;j

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

四阶龙格库塔法解微分方程 一、四阶龙格库塔法解一阶微分方程 ①一阶微分方程：cos ，初始值y(0)=0，求 y t 解区间[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; %%%%% 初始值

控制系统数字仿真 四阶龙格库塔法

控制系统数字仿真 1.实验目的 1.掌握利用四阶龙格-库塔（Runge-Kutta）法进行控制系统数字仿真的方 法。 2.学习分析高阶系统动态性能的方法。 3.学习系统参数改变对系统性能的影响。 二、实验内容 已知系统结构如下图 若输入为单位阶跃函数，计算当超调量分别为5%，25%，和50%时K的取值（用主导极点方法估算），并根据确定的K值在计算机上进行数字仿真。 三、实验过程 1.计算K值 二阶系统单位阶跃响应的超调量 %100% =? 1.当σ%=5%时

解得 ζ=0.690 设主导极点 =ζa + a=0.69a+j0.72a 代入D （s ）= 32 1025s s s K +++=0中， 32(0.690.72)10(0.690.72)25(0.690.72)0 a j a a j a a j a K ++++++=解得K=31.3，a=-2.10 即1,2 1.45 1.52s j =-± 2. 当σ%=25%时 解得 ζ=0.403 设主导极点 =ζa + a=0.403a+j0.915a 代入D （s ）= 321025s s s K +++=0中， 32(0.4030.915)10(0.4030.915)25(0.4030.915)0 a j a a j a a j a K ++++++=解得K=59.5，a=-2.75 即1,2 1.11 2.53s j =-± 3. 当σ%=50%时 解得 ζ=0.215 设主导极点 =ζa + a=0.215a+j0.977a 代入D （s ）= 321025s s s K +++=0中， 32(0.2150.977)10(0.2150.977)25(0.2150.977)0 a j a a j a a j a K ++++++=解得K=103，a=-3.48

龙格库塔法求微分方程2

《MATLAB 程序设计实践》课程考核 一、编程实现“四阶龙格－库塔（R-K ）方法求常微分方程”，并举一 例应用之。 【实例】采用龙格-库塔法求微分方程： ?? ?==+-=0 , 0)(1 '00 x x y y y 1、算法说明： 在龙格-库塔法中，四阶龙格-库塔法的局部截断误差约为o(h5)，被广泛应用于解微分方程的初值问题。其算法公式为： )22(6 3211k k k h y y n n +++=+ 其中： ?????????++=++=++ ==) ,() 21 ,21()21 ,21() ,(34 23121hk 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<=0 PointNum=100; end if nargin<3 y0=0; end y(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) end 3.2、实例求解源程序： %运行四阶R-K法

蒙特卡洛求定积及龙哥库塔解微分方程

蒙特卡洛法求定积分： 整体思路，蒙特卡洛求定积分的主要思想就是通过在取值范围内大量随机数的随机选取对函数进行求值，进而除以所取次数并乘以其区间，即为所积值。 Step 1： 在执行程序前，打开matlab，执行以下操作打开File—>New—>Function,并点击Function，弹出Editor窗口，将其中内容修改为 function [ y ] = f( x ) y=cos(x)+2.0; end （如图所示）。 Step 2： 在Editor窗口中点击File—>Save As，保存至所建立的文件夹内，保存名称必须为f.m。 Step 3： 回到Matlab程序中，将Current Folder改为与刚刚Function函数定义的保存路径一致的路径。Step 4： 在Command Window里输入以下源程序。 源程序： >> n=10^6; %用来表示精度，表示需要执行次数。 >> y=0; %初始化y=0，因为f(x)=cos(x)+2.0，下面出现y=y+f(x),故尔y用来统计计算出的所有f(x)值的和。 >> count=1; %从1开始计数，显然要执行n次。 >> while count<=n %当执行次数少于n次时，继续执行 x=unifrnd(0,4); %在matlab中，unifrnd是在某个区间内均匀选取实数（可为小数或整 数），表示x取0到4的任意数。 y=y+f(x); %求出到目前为止执行出的所有f(x)值的和，并用y表示 count=count+1; %count+1表示已执行次数再加一，将来通过与n的比对，判断是否执行下一次 end %如果执行次数已达n次，那么结束循环 z=y/n*4 %用y除以执行次数n,那么取得平均值,并用它乘以区间长度4,即可得到z 的近似值 z=7.2430 由于matlab中不能使用θ字符，故用z代替，即可得到θ=7.2395。 在执行过程中，发现每一次执行结果都会不一样，这是因为是随机选取，所以不同次数的计算结果会有偏差，但由于执行次数很多，从概率的角度来讲都是与真实值相近似的值。

四阶龙格库塔法的编程(赵)

例题一 四阶龙格-库塔法的具体形式为： 1.1程序: /*e.g: y'=t-y,t∈[0,1] /*y(0)=0 /*使用经典四阶龙格-库塔算法进行高精度求解 /* y(i＋1)＝yi＋( K1+ 2*K2 +2*K3+ K4)/6 /* K1=h*f(ti,yi) /* K2=h*f(ti+h/2,yi+K1*h/2) /* K3=h*f(ti+h/2,yi+K2*h/2) /* K4=h*f(ti+h,yi+K3*h) */ #include #include #define N 20 float f(float d,float p) //要求解的微分方程的右部的函数e.g: t-y { float derivative; derivative=d-p; return(derivative); } void main() { float t0; //范围上限

float t; //范围下限 float h; //步长 float nn; //计算出的点的个数,即迭代步数 int n; //对nn取整 float k1,k2,k3,k4; float y[N]; //用于存放计算出的常微分方程数值解 float yy; //精确值 float error;//误差 int i=0,j; //以下为函数的接口 printf("input t0:"); scanf("%f",&t0); printf("input t:"); scanf("%f",&t); printf("input y[0]:"); scanf("%f",&y[0]); printf("input h:"); scanf("%f",&h); // 以下为核心程序 nn=(t-t0)/h; printf("nn=%f\n",nn); n=(int)nn; printf("n=%d\n",n); for(j=0;j<=n;j++) { yy=t0-1+exp(-t0); //解析解表达式 error=y[i]-yy; //误差计算 printf("y[%f]=%f yy=%f error=%f\n",t0,y[i],yy,error);//结果输出k1=h*f(t0,y[i]); //求K1 k2=h*f((t0+h/2),(y[i]+k1*h/2)); //求K2 k3=h*f((t0+h/2),(y[i]+k2*h/2)); //求K3

龙格-库塔法求微分方程2

《MATLAB 程序设计实践》课程考核 一、编程实现“四阶龙格－库塔（R-K ）方法求常微分方程”，并举一 例应用之。 【实例】采用龙格-库塔法求微分方程： ? ? ?==+-=0 , 0)(1 '00x x y y y 1、算法说明： 在龙格-库塔法中，四阶龙格-库塔法的局部截断误差约为o(h5)，被广泛应用于解微分 方程的初值问题。其算法公式为： )22(6 3211k k k h y y n n +++ =+ 其中： ?????????++=+ + =++==) ,() 21 ,21()21 ,21(),(34 23121hk 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<=0 PointNum=100; end if nargin<3 y0=0; end y(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) end 3.2、实例求解源程序： %运行四阶R-K法

四阶龙格-库塔法解微分方程(C++)

一.作业： 用四阶龙格—库塔法求下列方程： ()()()1010100y x y x y '=-≤≤???=??步长0.1h =。并与解析解2 51x y e -=-比较。 二.程序 #include #include #include using namespace std; double f(double x,double y){ return 10*x*(1-y); } int main(){ int n,i; double h,k1,k2,k3,k4; cout<<"Please input the number of intervals:"; cin>>n; double *x=new double [n+1]; double *y=new double [n+1]; double *y1=new double [n+1]; h=1.0/n; y[0]=0; for(i=0;i<=n;i++){ x[i]=i*h; k1=f(x[i],y[i]); k2=f(x[i]+h/2,y[i]+k1*h/2); k3=f(x[i]+h/2,y[i]+k2*h/2); k4=f(x[i]+h,y[i]+k3*h); y[i+1]=y[i]+(k1+2*k2+2*k3+k4)*h/6; y1[i]=1-exp(-5*pow(x[i],2)); cout<

左侧为数值解,右侧为解析解.

四阶龙格库塔法原理C代码

/** ***四阶Runge-Kutta法*** 经典格式： y(n+1) = y(n) + h/6 ( K1 + 2*K2 + 2*K3 + K4 ) K1 = f( x(n) , y(n) ) K2 = f( x(n+1/2) , y(n) + h/2*K1 ) K3 = f( x(n+1/2) , y(n) + h/2*K2 ) K4 = f( x(n+1) , y(n) + h*K3 ) Runge-Kutta法是基于泰勒展开方法，因而需要所求解具有较好的光滑性。 属性：差分方法 《数值分析简明教程》-2 Editon -高等教育出版社- page 105 算法流程图代码维护：2005.6.14 DragonLord **/ #include #include #include /* 举例方程： y'= y - 2*x / y ( 0>x0>>y0>>h>>N) { int n=0;

for(;n

求解常微分方程组初值问题的龙格库塔法分析及其C代码

求解常微分方程组初值问题的 龙格库塔法分析及其C 代码 1、概 述 由高等数学的知识可知，一些特殊类型的常微分方程（组）能够求出给定初始值的解析解，而在科学与工程问题中遇到的常微分方程（组）往往是极其复杂的，要想求得其给定初始值的解析解就变得极其困难，甚至是得不到解析解。尽管如此，在研究实际问题时，往往只需要获得若干点上的近似值就行了，这就说明研究常微分方程（组）的数值解法是很有必要的。求解常微分方程（组）的数值解法有多种，比如欧拉法、龙格库塔法、线性多步法等等，其中龙格库塔法是这几种方法中比较常用的，因为其计算精度较高且具有自启动的特点。 对于用龙格库塔法求解单个常微分方程和求解常微分方程组的思路基本相似（注意一点一个微分方程组是常微分方程组即表明微分方程中的各阶导数都是对同一个变量求导，例如可以把各个量对时间求导得到一个常微分方程组，如果一个微分方程组中的有对不同变量的导数那么这个方程组就成了偏微分方程组），都是根据泰勒展开得到其迭代计算形式，其基本思想都是按照一定原则求取当前点附近一些点的斜率，通过这些斜率的线性组合作为当前点处的斜率，进行递推求解。 2、数学模型 2.1 常微分方程初值问题的数学模型 000 (,) ()()dy f x y x x dx y x y ?=>???=? （1） 式中(,) f x y 为,x y 的已知函数，0y 为给定的初始值。 常微分方程的数值解法的任务就是要求出函数值y 在微分自变量x 取如下序列： 012n x x x x <<< 时的值() (1,2,3,)i y x i n = 的近似值(1,2,3,)i y i n = ，一般情况下都采取等距结点的方式，即：0 (1,2,)i x x ih i n =+= ，其中h 为相邻两结点的距离，称为步长。 2.2 常微分方程组初值问题的数学模型

标准四阶龙格——库塔法

实验名：常微分方程数值解法 实习目的： （1） 通过实习进一步掌握标准四阶龙格——库塔法的基本思想； （2） 通过对标准四阶龙格——库塔法的调试练习，进一步体会其特点； （3） 通过实习进一步掌握标准四阶龙格——库塔法的计算步骤，并能灵活应用； （4） 通过上机调试运行，逐步培养解决实际问题的编程能力。 实习要求： （1） 熟悉Turbo C 的编译环境； （2） 实习前复习标准四阶龙格——库塔法的基本思想和过程； （3） 实习前复习标准四阶龙格——库塔法的计算步骤。 实习设备： （1） 硬件设备：单机或网络环境下的微型计算机一台; （2） 软件设备：DOS3.3以上操作系统，Turbo C2.0编译器。 实习内容： 标准四阶龙格——库塔法： （1）使用标准四阶龙格——库塔法求解初值问题 的数值求解。 （2）要求： 请写出程序的运行结果： 程序代码： #include "stdio.h" #include "conio.h" float func(float x,float y) { return(2*x*y); } float runge_kutta(float x0,float xn,float y0,int n) { float x,y,y1,y2,h,xh; float d1,d2,d3,d4; int i; x=x0; y=y0; h=(xn-x0)/n; for(i=1;i<=n;i++) 1(0)y 1x 0 2xy y =≤≤='

{ xh=x+h/2; d1=func(x,y); d2=func(xh,y+h*d1/2.0); d3=func(xh,y+h*d2/2.0); d4=func(xh,y+h*d3); y=y+h*(d1+2*d2+2*d3+d4)/6.0; x=x0+i*h; } return(y); } void main() { float x0,xn,y0,e; int n; printf("\ninput n:\n"); scanf("%d",&n); printf("input x0,xn:\n"); scanf("%f%f",&x0,&xn); printf("input y0:\n"); scanf("%f",&y0); e=runge_kutta(x0,xn,y0,n); printf("y(%f)=%6.6f",y0,e); } 运行结果: （3）思考题： 标准四阶龙格——库塔法的基本思想是什么? 龙格和库塔提出了一种间接地运用Taylor公式的方法，即利用y(x)在若干个待定点上的函数值和导数值做出线性组合式，选取适当系数使这个组合式进Taylor展开后与y(xi+1)的Taylor 展开式有较多的项达到一致，从而得出较高阶的数值公式，这就是龙格—库塔法的基本思想。