龙格库塔法求微分方程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法

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σ），并求出该应力状态下任意方位角α的斜截面上的应力σ、τ。

2、程序流程图：

σ

x

σy

σy

σx

σα

τ

yx

τyx

τxy

τxy

τ

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('应力圆');