中南大学
MATLAB 程序设计实践
材料科学与工程学院 2013年3月26日
一、编程实现“四阶龙格－库塔（R-K ）方法求常微分方程”，并举一
例应用之。
【实例】采用龙格-库塔法求微分方程：
⎧y ' =-y +1
⎨
⎩y (x 0) =0 , x 0=0
1、算法说明：
在龙格-库塔法中，四阶龙格-库塔法的局部截断误差约为o(h5)，被广泛应用于解微分方程的初值问题。其算法公式为：
h
y n +1=y n +(k 1+2k 2+2k 3)
6
其中：
⎧k 1=f (x n , y n ) ⎪
⎪k 2=f (x n +1h , y n +1hk 1) ⎪22 ⎨
⎪k 3=f (x n +1h , y n +1hk 2) ⎪22⎪k =f (x +h , y +hk )
n n 3⎩4
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
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 圆。
τyx
σy
σ
τ
τxy
σx
α
τxy
σx
τyx
σy
1、程序说明：
利用平面应力状态下斜截面应力的一般公式，画出任意平面应力状态下的应力圆
（Moore 圆），求出相应平面应力状态下的主应力（σ1、σ3），并求出该应力状态下任意方位角α的斜截面上的应力σ、τ。
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
MPa, σy =30 MPa, τxy =-20 MPa 为例） 4、程序运行结果：（以σx =100
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) (x -2)+(y -3+2x )=5 ; (2) 2(x -3)+(y /3)=4
2
2
2
2
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]; intersection_x_y=double(middle)
4、程序运行结果：
%画第一个椭圆 %画第二个椭圆 %显示网格 %输入两个椭圆方程 %联立两个椭圆方程求解交点 %合并x,y 两个矩阵 %将符号解转换成数值解
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
中南大学
MATLAB 程序设计实践
材料科学与工程学院 2013年3月26日
一、编程实现“四阶龙格－库塔（R-K ）方法求常微分方程”，并举一
例应用之。
【实例】采用龙格-库塔法求微分方程：
⎧y ' =-y +1
⎨
⎩y (x 0) =0 , x 0=0
1、算法说明：
在龙格-库塔法中，四阶龙格-库塔法的局部截断误差约为o(h5)，被广泛应用于解微分方程的初值问题。其算法公式为：
h
y n +1=y n +(k 1+2k 2+2k 3)
6
其中：
⎧k 1=f (x n , y n ) ⎪
⎪k 2=f (x n +1h , y n +1hk 1) ⎪22 ⎨
⎪k 3=f (x n +1h , y n +1hk 2) ⎪22⎪k =f (x +h , y +hk )
n n 3⎩4
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
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 圆。
τyx
σy
σ
τ
τxy
σx
α
τxy
σx
τyx
σy
1、程序说明：
利用平面应力状态下斜截面应力的一般公式，画出任意平面应力状态下的应力圆
（Moore 圆），求出相应平面应力状态下的主应力（σ1、σ3），并求出该应力状态下任意方位角α的斜截面上的应力σ、τ。
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
MPa, σy =30 MPa, τxy =-20 MPa 为例） 4、程序运行结果：（以σx =100
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) (x -2)+(y -3+2x )=5 ; (2) 2(x -3)+(y /3)=4
2
2
2
2
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]; intersection_x_y=double(middle)
4、程序运行结果：
%画第一个椭圆 %画第二个椭圆 %显示网格 %输入两个椭圆方程 %联立两个椭圆方程求解交点 %合并x,y 两个矩阵 %将符号解转换成数值解
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