四阶龙格库塔RK方法求常微分方程Word格式文档下载.docx
《四阶龙格库塔RK方法求常微分方程Word格式文档下载.docx》由会员分享，可在线阅读，更多相关《四阶龙格库塔RK方法求常微分方程Word格式文档下载.docx（13页珍藏版）》请在冰豆网上搜索。
%自变量在[x0,xt]上取的点数：
PointNum
%varargin为可输入项.可传适当参数给函数f（x,y）
%x：
所取的点的x值
%y：
对应点上的函数值
ifnargin<
4|PointNum<
=0
PointNum=100;
end
3
y0=0;
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）
3.2、实例求解源程序：
%运行四阶R-K法
clear,clc%清除内存中的变量
x0=0;
xt=2;
Num=100;
h=（xt-x0）/（Num-1）;
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'
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,'
%画出应力圆.标出该应力状态下各应力参数
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='
%控制坐标轴颜色为黑色
line（[b
（1）,b
（2）],[0,0],ps_array）;
%调整坐标轴
line（[0,0],[b（3）,b（4）],ps_array）;
%画出x坐标轴
plot（Sa,0,'
）%标出圆心
text（Sa,0,'
O'
plot（Smax,0,'
Smin,0,'
Sa,Tmax,'
Sa,Tmin,'
）%标出最大、最小拉应力、切应力点
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：
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
（以
为例）
Sigma_x（MPa）=100
Sigma_y（MPa）=30
Tau_xy（MPa）=-20
1
30
sigma=
99.8205
tau=
20.3109
Sigma_Max=
105.3087
Sigma_Min=
24.6970
Tau_Max=
40.3109
Tau_Min=
-40.2963
Sigma1=
Sigma3=
ans=
alpha_0=
14.8724
（二）实验5（椭圆的交点）两个椭圆可能具有0~4个交点.求下列两个椭圆的所有交点坐标：
（1）
（2）
此题相当于求两一个二元二次方程组的解.故为便于清楚地显示出两椭圆的相对位置.用ezplot函数把两个椭圆画在同一个坐标系中.然后利用fsolve函数解方程组得到两椭圆的交点即方程组的解。
clc;
ezplot（'
（x-2）^2+（y-3+2*x）^2-5'
[-1,5,-8,8]）;
%画第一个椭圆
2*（x-3）^2+（y/3）^2-4'
%画第二个椭圆
gridon;
%显示网格
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,'
%联立两个椭圆方程求解交点
middle=[x,y];
%合并x,y两个矩阵
intersection_x_y=double（middle）%将符号解转换成数值解
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
升级会员 