四阶龙格库塔RK方法求常微分方程.docx
- 文档编号:26075508
- 上传时间:2023-06-17
- 格式:DOCX
- 页数:12
- 大小:173.33KB
四阶龙格库塔RK方法求常微分方程.docx
《四阶龙格库塔RK方法求常微分方程.docx》由会员分享,可在线阅读,更多相关《四阶龙格库塔RK方法求常微分方程.docx(12页珍藏版)》请在冰豆网上搜索。
四阶龙格库塔RK方法求常微分方程
四阶龙格-库塔(R-K)方法求常微分方程
中南大学
MATLAB程序设计实践
材料科学与工程学院
2013年3月26日
一、编程实现“四阶龙格-库塔(R-K)方法求常微分方程”,并举一例应用之。
【实例】采用龙格-库塔法求微分方程:
1、算法说明:
在龙格-库塔法中,四阶龙格-库塔法的局部截断误差约为o(h5),被广泛应用于解微分方程的初值问题。
其算法公式为:
其中:
2、流程图:
2.1、四阶龙格-库塔(R-K)方法流程图:
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)
holdon
plot(x1,y1,'bo',xr,yr,'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,'.r');%画出应力圆,标出该应力状态下各应力参数
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='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坐标轴
holdon
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:
');
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
4、程序运行结果:
(以
为例)
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)
;
(2)
1、算法说明:
此题相当于求两一个二元二次方程组的解,故为便于清楚地显示出两椭圆的相对位置,用ezplot函数把两个椭圆画在同一个坐标系中,然后利用fsolve函数解方程组得到两椭圆的交点即方程组的解。
2、程序流程图:
3、程序代码:
clear;clc;
ezplot('(x-2)^2+(y-3+2*x)^2-5',[-1,5,-8,8]);%画第一个椭圆
holdon
ezplot('2*(x-3)^2+(y/3)^2-4',[-1,5,-8,8]);%画第二个椭圆
gridon;%显示网格
holdoff
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];%合并x,y两个矩阵
intersection_x_y=double(middle)%将符号解转换成数值解
4、程序运行结果:
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 方法 微分方程