中南大学MATLAB程序设计实践
材料科学与工程学院
2013年3月26日
一、编程实现“四阶龙格－库塔（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
2
3121hk 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法
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 on
plot(x1,y1,'bo',xr,yr,'r*')
4、程序运行结果：
MyRunge_Kutta_x =
0 0.5000 1.0000 1.5000 2.0000
MyRunge_Kutta_y =
0 0.3932 0.6318 0.7766 0.8645
二、变成解决以下科学计算问题：
（一）[例7-2-4] 材料力学复杂应力状态的分析——Moore 圆。
1、程序说明：
利用平面应力状态下斜截面应力的一般公式，画出任意平面应力状态下的应力圆
（Moore 圆），求出相应平面应力状态下的主应力（1σ、3σ），并求出该应力状态下任意方位角α的斜截面上的应力σ、τ。
σ
x
σy
σy
σx
σα
τ
yx
τyx
τxy
τxy
τ
2、程序流程图：
3、程序代码：
clear;clc;
Sx=input('Sigma_x(MPa)='); %输入该应力状态下的各应力值
Sy=input('Sigma_y(MPa)=');
Txy=input('Tau_xy(MPa)=');
a=linspace(0,pi,100); %等分半圆周角
Sa=(Sx+Sy)/2;
Sd=(Sx-Sy)/2;
Sigma=Sa+Sd*cos(2*a)-Txy*sin(2*a); %应力圆一般方程Tau=Sd*sin(2*a)+Txy*cos(2*a);
plot(Sigma,Tau,Sx,Txy,'.r',Sy,-Txy,'.r'); %画出应力圆，标出该应力状态下各应力参数
line([Sx,Sy],[Txy,-Txy]);
axis equal; %控制各坐标轴的分度使其相等——使应力圆变圆
title('应力圆');
xlabel('正应力（MPa）');
ylabel('剪应力（MPa）');
text(Sx,Txy,'A');
text(Sy,-Txy,'B');
Smax=max(Sigma);
Smin=min(Sigma);
Tmax=max(Tau);
Tmin=min(Tau);
b=axis; %提取坐标轴边界
ps_array.Color='k'; %控制坐标轴颜色为
黑色line([b(1),b(2)],[0,0],ps_array); %调整坐标轴
line([0,0],[b(3),b(4)],ps_array);
b=axis; %提取坐标轴边界line([b(1),b(2)],[0,0],ps_array); %画出x坐标轴
hold on
plot(Sa,0,'.r') %标出圆心
text(Sa,0,'O');
plot(Smax,0,'.r',Smin,0,'.r',Sa,Tmax,'.r',Sa,Tmin,'.r') %标出最大、最小拉
应力、切应力点text(Smax,0,'C');
text(Smin,0,'D');
text(Sa,Tmax,'E');
text(Sa,Tmin,'F');
%-----------此部分求某一斜截面上的应力---------------------------------------------
t=input('若需求某一截面上的应力，请输入1；若不求应力，请输入0：');
while t~=0
alpha=input('给出斜截面方向角：alpha=（角度）：');
sigma=Sa+Sd*cos(2*(alpha/180*pi))-Txy*sin(2*(alpha/180*pi))
tau=Sd*sin(2*(alpha/180*pi))+Txy*cos(2*(alpha/180*pi))
plot(sigma,tau,'or')
t=input('若还需求其他截面上的应力，请输入1；若要退出，请输入0：');
end
hold off
%----------此部分求该应力状态下的主应力--------------------------------------------
Sigma_Max=Smax
Sigma_Min=Smin
Tau_Max=Tmax
Tau_Min=Tmin
Sigma1=Smax %得出主应力 Sigma3=Smin
l=Sx-Sa;h=Txy;ratio=abs(h/l); %求主应力平面方向角 '主应力平面方向角（角度）：' alpha_0=atan(ratio)/2*180/pi
4、程序运行结果：（以MPa 20 MPa, 30 MPa, 100-===xy y x τσσ为例）
Sigma_x(MPa)=100 Sigma_y(MPa)=30 Tau_xy(MPa)=-20
若需求某一截面上的应力，请输入1；若不求应力，请输入0：1 给出斜截面方向角：alpha=（角度）：30 sigma = 99.8205
tau =
20.3109
若还需求其他截面上的应力，请输入1；若要退出，请输入0：0 Sigma_Max = 105.3087
Sigma_Min = 24.6970
Tau_Max = 40.3109
Tau_Min = -40.2963
Sigma1 = 105.3087
Sigma3 = 24.6970
ans =
主应力平面方向角（角度）：
alpha_0 = 14.8724
（二）实验5（椭圆的交点） 两个椭圆可能具有0 ~ 4个交点，求下列两个椭圆
的所有交点坐标：
(1) ()()52322
2
=+-+-x y x ; (2) ()()43/322
2
=+-y x
1、算法说明：
此题相当于求两一个二元二次方程组的解，故为便于清楚地显示出两椭圆的相对位置，用ezplot 函数把两个椭圆画在同一个坐标系中，然后利用fsolve 函数解方程组得到两椭圆的交点即方程组的解。
2、程序流程图：
3、程序代码：
clear; clc;
ezplot('(x-2)^2+(y-3+2*x)^2-5',[-1,5, -8,8]); %画第一个椭圆
hold on
ezplot('2*(x-3)^2+(y/3)^2-4',[-1,5, -8,8]); %画第二个椭圆
grid on; %显示网格
hold off
f1=sym('(x-2)^2+(y-3+2*x)^2=5'); %输入两个椭圆方程f2=sym('2*(x-3)^2+(y/3)^2=4');
[x,y]=solve(f1,f2,'x','y'); %联立两个椭圆方程求解交点middle=[x,y]; %合并x,y两个矩阵intersection_x_y=double(middle) %将符号解转换成数值解
4、程序运行结果：
intersection_x_y =
4.0287 - 0.0000i -4.1171 + 0.0000i 3.4829 + 0.0000i -
5.6394 + 0.0000i 1.7362 - 0.0000i -2.6929 + 0.0000i 1.6581 + 0.0000i 1.8936 - 0.0000i
《数值分析》课程实验报告 一、实验目的 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) %将得到结果 四、实验小结 通过实验结果分析可知，龙格库塔方法求解常微分方程能获得比较好的数值解，龙格库塔方法的数值解的精度较高，方法比较简便易懂。 姓名：樊元君学号：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阶龙格库塔法解微分方程的程序 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实现~ function ff=rk(yy,x0,y0,h,a,b)%yy为y的导函数,x0,y0,为初值,h为步长,a,b为区间 c=(b-a)/h+1;i1=1; %c为迭代步数；i1为迭代步数累加值 y=y0;z=zeros(c,6); %z生成c行，5列的零矩阵存放结果； %每行存放c次迭代结果，每列分别存放k1~k4及y的结果 for x=a:h:b if i1<=c k1=feval(yy,x,y); k2=feval(yy,x+h/2,y+(h*k1)/2); k3=feval(yy,x+h/2,y+(h*k2)/2); k4=feval(yy,x+h,y+h*k3); y=y+(h/6)*(k1+2*k2+2*k3+k4); z(i1,1)=x;z(i1,2)=k1;z(i1,3)=k2;z(i1,4)=k3;z(i1,5)=k4;z(i1,6)=y; i1=i1+1; end end fprintf(‘结果矩阵，第一列为x(n),第二列~第五列为k1~k4，第六列为y（n+1）的结果') z %在命令框输入下列语句 %yy=inline('x+y'); %>> rk(yy,0,1,0.2,0,1) %将得到结果 %结果矩阵，第一列为x(n),第二列~第五列为k1~k4第六列为y（n+1）的结果 %z = % 0 1.0000 1.2000 1.2200 1.4440 1.2428 % 0.2000 1.4428 1.6871 1.7115 1.9851 1.5836 % 0.4000 1.9836 2.2820 2.3118 2.6460 2.0442 % 0.6000 2.6442 3.0086 3.0451 3.4532 2.6510 % 0.8000 3.4510 3.8961 3.9407 4.4392 3.4365 % 1.0000 4.4365 4.9802 5.0345 5.6434 4.4401 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 函数功能编辑本段回目录 ode是专门用于解微分方程的功能函数，他有ode23,ode45,ode23s等等，采用的是Runge-Kutta算法。ode45表示采用四阶，五阶runge-kutta单步算法,截断误差为(Δx)3。解决的是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 龙格-库塔方法是一种经典方法，具有很高的精度，它间接的利用了泰勒级数展开，避免了高阶偏导数的计算。此处以最为经典的四级四阶龙格-库塔方法为例，计算格式如下 ()()()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,’*’); 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 Matlab中龙格-库塔(Runge-Kutta)方法原理及实现 龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法。由于此算法精度高，采取措施对误差进行抑制，所以其实现原理也较复杂。该算法是构建在数学支持的基础之上的。龙格库塔方法的理论基础来源于泰勒公式和使用斜率近似表达微分，它在积分区间多预计算出几个点的斜率，然后进行加权平均，用做下一点的依据，从而构造出了精度更高的数值积分计算方法。如果预先求两个点的斜率就是二阶龙格库塔法，如果预先取四个点就是四阶龙格库塔法。一阶常微分方程可以写作：y'=f(x,y),使用差分概念。 (Yn+1-Yn)/h= f(Xn,Yn)推出（近似等于，极限为Yn'） Yn+1=Yn+h*f(Xn,Yn) 另外根据微分中值定理，存在0 所以，为了更好更准确地把握时间关系，应自己在理解龙格库塔原理的基础上，编写定步长的龙格库塔函数，经过学习其原理，已经完成了一维的龙格库塔函数。 仔细思考之后，发现其实如果是需要解多个微分方程组，可以想象成多个微分方程并行进行求解，时间，步长都是共同的，首先把预定的初始值给每个微分方程的第一步，然后每走一步，对多个微分方程共同求解。想通之后发现，整个过程其实很直观，只是不停的逼近计算罢了。编写的定步长的龙格库塔计算函数： 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:n x(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; %按照龙格库塔方法进行数值求解 /*************************************************************************** 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 四阶龙格库塔法解微分方程 一、四阶龙格库塔法解一阶微分方程 ①一阶微分方程：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. 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:n x(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; %按照龙格库塔方法进行数值求解 end 2.另外再新建一个.,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，间隔为1 subplot(122) plot(T1,F1)%自编的龙格库塔函数效果 title('自编的龙格库塔函数') grid 运行步骤3文件即可得到结果，F1为估测值 或者可以调用matlab自带函数ode45求解 方法如下： 控制系统数字仿真 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 《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。 在执行过程中，发现每一次执行结果都会不一样，这是因为是随机选取，所以不同次数的计算结果会有偏差，但由于执行次数很多，从概率的角度来讲都是与真实值相近似的值。 龙格—库塔方法求解微分方程初值问题 （数学1201+41262022+陈晓云） 初值问题： y x x -+=2dx dy ，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 f K h n h n 232,2 ()K h y h x f K n n 34,++= ()K K K K y y h n 4 3211n 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)-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*(k1+k2*2+2*k3+k4)/6; end x=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 long y2=x.^2-x+1 R4=y2-y1 [x',y1',y2',R4'] y2=dsolve('Dy=x^2+x-y','y(0)=1','x') plot(x,y1,'b*-') hold on y3=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.840000000000000 0.900000000000000 0.910001299159352 0.910000000000000 1.000000000000000 1.000001383861433 1.000000000000000MATLAB改进欧拉法与四阶龙格-库塔求解一阶常微分方程
matlab编的4阶龙格库塔法解微分方程的程序
龙格库塔方法matlab实现
matlab 四阶龙格-库塔法求微分方程
经典四阶龙格库塔法解一阶微分方程组
Matlab中龙格-库塔(Runge-Kutta)方法原理及实现
MATLAB龙格-库塔方法解微分方程
龙格库塔算法解微分方程组-c语言
Matlab中龙格-库塔(Runge-Kutta)方法原理及实现
龙格库塔算法解微分方程组 c语言
四阶龙格库塔法解微分方程
龙格-库塔法MATLAB
控制系统数字仿真 四阶龙格库塔法
龙格库塔法求微分方程2
蒙特卡洛求定积及龙哥库塔解微分方程
龙格库塔法求微方程matlab
相关主题
文本预览