四阶龙格库塔RK方法求常微分方程Word格式.docx
- 文档编号:22480953
- 上传时间:2023-02-04
- 格式:DOCX
- 页数:13
- 大小:163.50KB
四阶龙格库塔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
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 四阶龙格库塔 RK 方法 微分方程