四阶龙格-库塔(R-K)方法求常微分方程

中南大学

四阶龙格-库塔(R-K)方法求常微分方程

MATLAB 程序设计实践

材料科学与工程学院 2012年4月10日

一、编程实现“四阶龙格-库塔(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 )方法流程图:

四阶龙格-库塔(R-K)方法求常微分方程

2.2、实例求解流程图:

四阶龙格-库塔(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<=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

四阶龙格-库塔(R-K)方法求常微分方程

二、变成解决以下科学计算问题:

(一)[例7-2-4] 材料力学复杂应力状态的分析——Moore 圆。

1、程序说明:

利用平面应力状态下斜截面应力的一般公式,画出任意平面应力状态下的应力圆

(Moore 圆),求出相应平面应力状态下的主应力(1σ、3σ),并求出该应力状态下任意方位角α的斜截面上的应力σ、τ。

σ

x

σy

σy

σx

σα

τ

yx

τyx

τxy

τxy

τ

2、程序流程图:

四阶龙格-库塔(R-K)方法求常微分方程

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

四阶龙格-库塔(R-K)方法求常微分方程

(二)实验5(椭圆的交点) 两个椭圆可能具有0 ~ 4个交点,求下列两个椭圆

的所有交点坐标:

(1) ()()52322

2

=+-+-x y x ; (2) ()()43/322

2

=+-y x

1、算法说明:

此题相当于求两一个二元二次方程组的解,故为便于清楚地显示出两椭圆的相对位置,用ezplot 函数把两个椭圆画在同一个坐标系中,然后利用fsolve 函数解方程组得到两椭圆的交点即方程组的解。

2、程序流程图:

四阶龙格-库塔(R-K)方法求常微分方程

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、程序运行结果:

四阶龙格-库塔(R-K)方法求常微分方程

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

相关推荐
相关主题
热门推荐