四阶龙格库塔RK方法求常微分方程.docx
《四阶龙格库塔RK方法求常微分方程.docx》由会员分享，可在线阅读，更多相关《四阶龙格库塔RK方法求常微分方程.docx（12页珍藏版）》请在冰豆网上搜索。
四阶龙格库塔RK方法求常微分方程
中南大学
MATLAB程序设计实践
材料科学与工程学院
2013年3月26日
一、编程实现“四阶龙格－库塔（R-K）方法求常微分方程”，并举一例应用之。
【实例】采用龙格-库塔法求微分方程：
1、算法说明：
在龙格-库塔法中，四阶龙格-库塔法的局部截断误差约为o（h5），被广泛应用于解微分方程的初值问题。
其算法公式为：
其中：
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：
对应点上的函数值
ifnargin<4|PointNum<=0
PointNum=100;
end
ifnargin<3
y0=0;
end
y（1,:
）=y0（:
）';%初值存为行向量形式
h=（xt-x0）/（PointNum-1）;%计算步长
x=x0+[0:
（PointNum-1）]'*h;%得x向量值
fork=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）
holdon
plot（x1,y1,'bo',xr,yr,'r*'）
4、程序运行结果：
MyRunge_Kutta_x=
00.50001.00001.50002.0000
MyRunge_Kutta_y=
00.39320.63180.77660.8645
二、变成解决以下科学计算问题：
（一）[例7-2-4]材料力学复杂应力状态的分析——Moore圆。
1、程序说明：
利用平面应力状态下斜截面应力的一般公式，画出任意平面应力状态下的应力圆（Moore圆），求出相应平面应力状态下的主应力（
、
），并求出该应力状态下任意方位角
的斜截面上的应力
、
。
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]）;
axisequal;%控制各坐标轴的分度使其相等——使应力圆变圆
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坐标轴
holdon
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：
'）;
whilet~=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
holdoff
%----------此部分求该应力状态下的主应力--------------------------------------------
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、程序运行结果：
（以
为例）
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）
;
（2）
1、算法说明：
此题相当于求两一个二元二次方程组的解，故为便于清楚地显示出两椭圆的相对位置，用ezplot函数把两个椭圆画在同一个坐标系中，然后利用fsolve函数解方程组得到两椭圆的交点即方程组的解。
2、程序流程图：
3、程序代码：
clear;clc;
ezplot（'（x-2）^2+（y-3+2*x）^2-5',[-1,5,-8,8]）;%画第一个椭圆
holdon
ezplot（'2*（x-3）^2+（y/3）^2-4',[-1,5,-8,8]）;%画第二个椭圆
gridon;%显示网格
holdoff
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.0000i1.8936-0.0000i
升级会员 