四阶龙格库塔RK方法求常微分方程文档格式.docx
《四阶龙格库塔RK方法求常微分方程文档格式.docx》由会员分享，可在线阅读，更多相关《四阶龙格库塔RK方法求常微分方程文档格式.docx（11页珍藏版）》请在冰豆网上搜索。
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：
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
（以
为例）
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
升级会员 