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

1、四阶龙格库塔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=inlin

2、e(-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=xrMyRunge_Kutta_y=yrplot(x,yt,k,x1,y1,b-,xr,yr,r-)legend(真值,ode23,Rung-Kutta法解,2)hold onplot(x1,y1,bo,xr,yr,r*)4、程序运行结果：MyRunge_Kutta_x = 0 0.5000 1.0000 1.

3、5000 2.0000MyRunge_Kutta_y = 0 0.3932 0.6318 0.7766 0.8645二、变成解决以下科学计算问题：（一）例7-2-4 材料力学复杂应力状态的分析Moore圆。1、程序说明：利用平面应力状态下斜截面应力的一般公式，画出任意平面应力状态下的应力圆（Moore圆），求出相应平面应力状态下的主应力（、），并求出该应力状态下任意方位角的斜截面上的应力、。2、程序流程图：3、程序代码：clear;clc;Sx=input(Sigma_x(MPa)=); %输入该应力状态下的各应力值Sy=input(Sigma_y(MPa)=);Txy=input(Tau_x

4、y(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);axis equal; %控制各坐标轴的分度使其相等使应力圆变圆title(应力圆);xlabel(正应力（MPa）);ylabel(剪应力（MPa）);text(Sx,

5、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坐标轴hold onplot(Sa,0,.r) %标出圆心text(Sa,0,O);plot(Smax,0,.r

6、,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：);while t=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

7、*cos(2*(alpha/180*pi) plot(sigma,tau,or) t=input(若还需求其他截面上的应力，请输入1；若要退出，请输入0：);endhold off%-此部分求该应力状态下的主应力-Sigma_Max=SmaxSigma_Min=SminTau_Max=TmaxTau_Min=TminSigma1=Smax %得出主应力Sigma3=Sminl=Sx-Sa;h=Txy;ratio=abs(h/l); %求主应力平面方向角主应力平面方向角（角度）：alpha_0=atan(ratio)/2*180/pi4、程序运行结果：（以为例）Sigma_x(MPa)=100S

8、igma_y(MPa)=30Tau_xy(MPa)=-20若需求某一截面上的应力，请输入1；若不求应力，请输入0：1给出斜截面方向角：alpha=（角度）：30sigma = 99.8205tau = 20.3109若还需求其他截面上的应力，请输入1；若要退出，请输入0：0Sigma_Max = 105.3087Sigma_Min = 24.6970Tau_Max = 40.3109Tau_Min = -40.2963Sigma1 = 105.3087Sigma3 = 24.6970ans =主应力平面方向角（角度）：alpha_0 = 14.8724（二）实验5（椭圆的交点） 两个椭圆可能具

9、有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); %画第一个椭圆hold onezplot(2*(x-3)2+(y/3)2-4,-1,5, -8,8); %画第二个椭圆grid on; %显示网格hold offf1=sym(x-2)2+(y-3+2*x

10、)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.0000i 1.8936 - 0.0000i