四阶龙格库塔RK方法求常微分方程Word格式文档下载.docx

1、%自变量在x0,xt上取的点数：PointNum%varargin为可输入项.可传适当参数给函数f（x,y）%x：所取的点的x值%y：对应点上的函数值 if nargin4 | PointNum=0 PointNum=100;end3 y0=0;y（1,:）=y0（:）; %初值存为行向量形式h=（xt-x0）/（PointNum-1）; %计算步长 x=x0+0:（PointNum-1）*h; %得x向量值for k=1:（PointNum） %迭代计算 f1=h*feval（fun,x（k）,y（k,:）,varargin:）; f1=f1（: %得公式k1 f2=h*feval（fun,

2、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

3、-exp（-a*x）; %真值解fun=inline（-y+1,xy）; %用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（真值ode23Rung-Kutta法解,2）hold onplot（x1,y1,bor*4、程序运行结果：MyRunge_Kutta_x = 0 0.

4、5000 1.0000 1.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=in

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

6、,Atext（Sy,-Txy,BSmax=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,Oplot（Smax,0,Smin,0,Sa,Tmax,Sa,Tmin,） %标出最大、最小拉应力、切应力点text（Smax,0,Ctext（Smin,0,Dt

7、ext（Sa,Tmax,Etext（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*cos（2*（alpha/180*pi） plot（sigma,tau,or t=input（若还需求其他截面上的应力.请输入1；若要退出.请输入0：hold off%-此部分

8、求该应力状态下的主应力-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/pi（以为例）Sigma_x（MPa）=100Sigma_y（MPa）=30Tau_xy（MPa）=-20130sigma = 99.8205tau = 20.3109Sigma_Max = 105.3087Sigma_Min = 24.6970Tau_M

9、ax = 40.3109Tau_Min = -40.2963Sigma1 =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 %画第二个椭圆grid on; %显

10、示网格f1=sym（x-2）2+（y-3+2*x）2=5 %输入两个椭圆方程f2=sym（2*（x-3）2+（y/3）2=4x,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.0000i 1.8936 - 0.0000i