四阶龙格库塔RK方法求常微分方程

《四阶龙格库塔RK方法求常微分方程》由会员分享，可在线阅读，更多相关《四阶龙格库塔RK方法求常微分方程（8页珍藏版）》请在装配图网上搜索。

1、中南大学MATLAB程序设计实践材料科学与工程学院20#3月26日一、编程实现四阶龙格库塔R-K方法求常微分方程,并举一例应用之.实例采用龙格-库塔法求微分方程：1、算法说明：在龙格-库塔法中,四阶龙格-库塔法的局部截断误差约为o,被广泛应用于解微分方程的初值问题.其算法公式为：其中：2、流程图：2.1、四阶龙格库塔R-K方法流程图：输入待求微分方程、求解的自变量范围、初值以与求解范围内的取点数等.确定求解范围内的步长k = 取点数？否求解：求解并输出：是结束算法2.2、实例求解流程图：输入求解的自变量范围求出待求简单微分方程的真值解用MATLAB自带函数ode23求解待求微分方程结束用自编函

2、数四阶龙格库塔R-K方法求解待求微分方程开始3、源程序代码3.1、四阶龙格库塔R-K方法源程序：function x,y = MyRunge_Kutta%Runge-Kutta 方法解微分方程形为 y=fx,y%此程序可解高阶的微分方程.只要将其形式写为上述微分方程的向量形式%函数 f: fun%自变量的初值和终值：x0, xt%y0表示函数在x0处的值,输入初值为列向量形式%自变量在x0,xt上取的点数：PointNum%varargin为可输入项,可传适当参数给函数f%x：所取的点的x值%y：对应点上的函数值if nargin4 | PointNum=0 PointNum=100;endi

3、f nargin3 y0=0;endy=y0; %初值存为行向量形式h=/; %计算步长 x=x0+0:*h; %得x向量值for k=1:%迭代计算f1=h*fevalfun,x,y,varargin:; f1=f1; %得公式k1 f2=h*fevalfun,x+h/2,y+f1/2,varargin:; f2=f2; %得公式k2 f3=h*fevalfun,x+h/2,y+f2/2,varargin:; f3=f3; %得公式k3 f4=h*fevalfun,x+h,y+f3,varargin:; f4=f4; %得公式k4 y=y+f1+2*+f4/6; %得yend3.2、实例求解

4、源程序：%运行四阶R-K法clear, clc %清除内存中的变量x0=0;xt=2;Num=100;h=/;x=x0+0:Num*h;a=1;yt=1-exp; %真值解fun=inline; %用inline构造函数fy0=0; %设定函数初值PointNum=5; %设定取点数x1,y1=ode23;xr,yr=MyRunge_Kutta;MyRunge_Kutta_x=xrMyRunge_Kutta_y=yrplotlegendhold onplot4、程序运行结果：MyRunge_Kutta_x = 0 0.5000 1.0000 1.5000 2.0000MyRunge_Kutta

5、_y = 0 0.3932 0.6318 0.7766 0.8645二、变成解决以下科学计算问题：一例7-2-4 材料力学复杂应力状态的分析Moore圆.1、程序说明：利用平面应力状态下斜截面应力的一般公式,画出任意平面应力状态下的应力圆Moore圆,求出相应平面应力状态下的主应力、,并求出该应力状态下任意方位角的斜截面上的应力、.2、程序流程图：开始输入待求应力状态的参数画出应力圆求某一方向角截面上的应力?输入方向角求出相应正应力、切应力是否得出该应力状态下的主应力求出主应力平面方向角结束3、程序代码：clear;clc;Sx=inputSigma_x=; %输入该应力状态下的各应力值Sy=

6、inputSigma_y=;Txy=inputTau_xy=;a=linspace; %等分半圆周角Sa=/2;Sd=/2;Sigma=Sa+Sd*cos-Txy*sin; %应力圆一般方程Tau=Sd*sin+Txy*cos;plot; %画出应力圆,标出该应力状态下各应力参数line;axis equal; %控制各坐标轴的分度使其相等使应力圆变圆title;xlabel;ylabel;text;text;Smax=max;Smin=min;Tmax=max;Tmin=min;b=axis; %提取坐标轴边界ps_array.Color=k; %控制坐标轴颜色为黑色lineb,b,0,0,

7、ps_array; %调整坐标轴line0,0,b,b,ps_array; b=axis; %提取坐标轴边界lineb,b,0,0,ps_array; %画出x坐标轴hold onplot %标出圆心text;plot %标出最大、最小拉应力、切应力点text;text;text;text;%-此部分求某一斜截面上的应力-t=input;while t=0 alpha=input; sigma=Sa+Sd*cos2*-Txy*sin2* tau=Sd*sin2*+Txy*cos2* plot t=input;endhold off%-此部分求该应力状态下的主应力-Sigma_Max=SmaxS

8、igma_Min=SminTau_Max=TmaxTau_Min=TminSigma1=Smax %得出主应力Sigma3=Sminl=Sx-Sa;h=Txy;ratio=abs; %求主应力平面方向角主应力平面方向角角度：alpha_0=atan/2*180/pi4、程序运行结果：以为例Sigma_x=100Sigma_y=30Tau_xy=-20若需求某一截面上的应力,请输入1；若不求应力,请输入0：1给出斜截面方向角：alpha=角度：30sigma = 99.8205tau = 20.3109若还需求其他截面上的应力,请输入1；若要退出,请输入0：0Sigma_Max = 105.30

9、87Sigma_Min = 24.6970Tau_Max = 40.3109Tau_Min = -40.2963Sigma1 = 105.3087Sigma3 = 24.6970ans =主应力平面方向角角度：alpha_0 = 14.8724二实验5椭圆的交点 两个椭圆可能具有0 4个交点,求下列两个椭圆的所有交点坐标： ; 1、算法说明：此题相当于求两一个二元二次方程组的解,故为便于清楚地显示出两椭圆的相对位置,用ezplot函数把两个椭圆画在同一个坐标系中,然后利用fsolve函数解方程组得到两椭圆的交点即方程组的解.2、程序流程图：开始依照所给方程画出两个椭圆用符号解法求两椭圆的交点将

10、符号解转换成数值解结束3、程序代码：clear; clc;ezplot2+2-5,-1,5, -8,8; %画第一个椭圆hold onezplot2*2+2-4,-1,5, -8,8; %画第二个椭圆grid on; %显示网格hold offf1=sym2+2=5; %输入两个椭圆方程f2=sym2*2+2=4;x,y=solve; %联立两个椭圆方程求解交点middle=x,y; %合并x,y两个矩阵intersection_x_y=double %将符号解转换成数值解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.0000i 1.8936 - 0.0000i8 / 8