四階龍格-庫塔(R-K)方法求常微分方程

1、實用標(biāo)準(zhǔn)文案中南大學(xué)MATLAB!序設(shè)計實踐材料科學(xué)與工程學(xué)院2013年3月26日一、編程實現(xiàn)“四階龍格庫塔（R-K）方法求常微分方程”，并舉一例應(yīng)用之?！緦嵗坎捎谬埜?庫塔法求微分方程：y' = -y +1y(x0) =0, X0 = 01、算法說明：在龍格-庫塔法中，四階龍格-庫塔法的局部截斷誤差約為 o（h5）,被廣泛應(yīng)用于解微分 方程的初值問題。其算法公式為：hyn 1 = yn (k1 2k2 2k3)6其中:k1 =f(Xn,yn)k2 = f (Xnk3 = f (Xnk4 = f(Xn1 .1 、2h, yn 2hk1)11-h, yn -hk2)2、流程圖：2.1、

2、四階龍格一庫塔（R-K）方法流程圖:輸入待求微分方程、求解的自變量范圍、 初值以及求解范圍內(nèi)的取點數(shù)等。kl =f (Xn, yn)11k2 =f(Xn j, yn +-hki) 22確定求解范圍內(nèi)的步長求解:h, yn hk3)11= f(Xn h, yn,. - hk2)22= f(Xn h,yn hk3)求解并輸出：yn十=yn +%.十羽十羽）2.2、實例求解流程圖:文檔3、源程序代碼3.1、 四階龍格一庫塔(R-K)方法源程序:function x,y = MyRunge_Kutta(fun,x0,xt,y0,PointNum,varargin)%Runge-Kutta方法解微分方程

3、形為y'(t)=f(x,y(x)%t程序可解高階的微分方程。只要將其形式寫為上述微分方程的向量形式%!數(shù) f(x,y): fun%變量的初值和終值：x0, xt%y成示函數(shù)在x0處的值，輸入初值為列向量形式%變量在x0,xt 上取的點數(shù)：PointNum%varargin為可輸入項，可傳適當(dāng)參數(shù)給函數(shù)f(x,y)%x所取的點的x值%y:對應(yīng)點上的函數(shù)值if nargin<4 | PointNum<=0 PointNum=100;endif nargin<3 y0=0;end%初值存為行向量形式%計算步長%4x向量值%迭代計算y(i,：)=y0(：);h=(xt-x0)

4、/(PointNum-1);x=x0+0:(PointNum-1)'*h;for k=1:(PointNum)f1=h*feval(fun,x(k),y(k,:),varargin:);f1=f1(:)'f2=h*feval(fun,x(k)+h/2,y(k,:)+f1/2,varargin:);f2=f2(:)'f3=h*feval(fun,x(k)+h/2,y(k,:)+f2/2,varargin:);f3=f3(:),;f4=h*feval(fun,x(k)+h,y(k,:)+f3,varargin:);f4=f4(:)'y(k+1,:)=y(k,:)+(

5、f1+2*(f2+f3)+f4)/6;end3.2、 實例求解源程序：%!行四階R-K法clear, clcx0=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');y0=0;PointNum=5;x1,y1=ode23(fun,0,2,0);xr,yr=MyRunge_Kutta(fun,x0,xt,y0,PointNum);MyRunge_Kutta_x=xr'%青除內(nèi)存中的變量%導(dǎo)公式k1%導(dǎo)

6、公式k2%導(dǎo)公式k3%導(dǎo)公式k4%4y(n+1)演值解%i inline構(gòu)造函數(shù) f(x,y)%設(shè)定函數(shù)初值%設(shè)定取點數(shù)MyRunge_Kutta_y=yr' plot(x,yt, 'k' ,x1,y1, 'b-' ,xr,yr, 'r-') legend(真值，ode23' , 'Rung-Kutta 法解',2) hold on plot(x1,y1, 'bo' ,xr,yr, 'r*')4、程序運行結(jié)果：MyRunge_Kutta_x =0 0,5000 1.0000 1.50

7、00 2.0000MyRunge_Kutta_y =0 0,3932 0.6318 0.7766 0.8645二、變成解決以下科學(xué)計算問題:（一）例 7-2-4材料力學(xué)復(fù)雜應(yīng)力狀態(tài)的分析Moore 圓。x1、程序說明：利用平面應(yīng)力狀態(tài)下斜截面應(yīng)力的一般公式， 畫出任意平面應(yīng)力狀態(tài)下的應(yīng)力圓 （Moore 圓），求出相應(yīng)平面應(yīng)力狀態(tài)下的主應(yīng)力（ 仃1、仃3 K并求出該應(yīng)力狀態(tài)下任意方位角 a的 斜截面上的應(yīng)力仃、丁。2、3、程序代碼:clear;clc;Sx=input( 'Sigma_x(MPa尸' 的各應(yīng)力值）；癡入該應(yīng)力狀態(tài)下Sy=input( 'Sigma_y(M

8、Pa尸'）；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' );%®出應(yīng)力圓，標(biāo)出該應(yīng)力狀態(tài)下各應(yīng)力參數(shù)line(Sx,Sy,Txy,-Txy);axis equal ;%控制各坐標(biāo)軸的分度使其相等一一使應(yīng)力

9、圓變圓title('應(yīng)力圓)；xlabel('正應(yīng)力(MPa ');ylabel('剪應(yīng)力(MPa ');text(Sx,Txy, 'A');text(Sy,-Txy, 'B');Smax=max(Sigma);Smin=min(Sigma);Tmax=max(Tau);%!取坐標(biāo)軸邊界%控制坐標(biāo)軸顏色為黑色%調(diào)整坐標(biāo)軸%提取坐標(biāo)軸邊界%畫出x坐標(biāo)軸%標(biāo)出圓心%標(biāo)出最大、最小拉應(yīng)力、切應(yīng)力點Tmin=min(Tau);b=axis;ps_array.Color= 'k'line(b(1),b(2),0,0,

10、ps_array);line(0,0,b(3),b(4),ps_array);b=axis;line(b(1),b(2),0,0,ps_array);hold onplot(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'

11、;);text(Sa,Tmin, 'F');% 此部分求某一斜截面上的應(yīng)力 t=input('若需求某一截面上的應(yīng)力，請輸入 1;若不求應(yīng)力，請輸入 0：'); while t=0alpha=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( '

12、;若還需求其他截面上的應(yīng)力，請輸入 1;若要退出，請輸入 0：');end hold off% 此部分求該應(yīng)力狀態(tài)下的主應(yīng)力 Sigma_Max=Smax Sigma_Min=Smin Tau_Max=Tmax Tau_Min=TminSigma1=SmaxSigma3=Sminl=Sx-Sa;h=Txy;ratio=abs(h/l);%得出主應(yīng)力%求主應(yīng)力平面方向角主應(yīng)力平面方向角(角度)：alpha_0=atan(ratio)/2*180/pi4、程序運行結(jié)果：(以仃* =100 MPa，仃v =30 MPa, txv =20 MPa為例) xyxySigma_x(MPa)=100

13、Sigma_y(MPa)=30Tau_xy(MPa)=-20若需求某一截面上的應(yīng)力，請輸入 1;若不求應(yīng)力，請輸入 0： 1給出斜截面方向角：alpha=(角度)：30sigma =99.8205tau =20.3109若還需求其他截面上的應(yīng)力，請輸入1;若要退出，請輸入 0: 0Sigma_Max =105.3087Sigma_Min =24.6970Tau_Max =40.3109Tau_Min =-40.2963Sigma1 =105.3087Sigma3 =24.6970 ans =主應(yīng)力平面方向角(角度)：alpha_0 =14.87240204060 BO 100正應(yīng)力(MPa)（

14、二）實驗5 （橢圓的交點）圓的所有交點坐標(biāo)：兩個橢圓可能具有0 4個交點，求下列兩個橢(1) x -2 2 y-3 2x 2 =5 ; (2)2x32y/3 2 = 41、算法說明：此題相當(dāng)于求兩一個二元二次方程組的解，故為便于清楚地顯示出兩橢圓的相對位置，用ezplot函數(shù)把兩個橢圓畫在同一個坐標(biāo)系中，然后利用fsolve函數(shù)解方程組得到兩橢圓的交點即方程組的解。2、程序流程圖：3、程序代碼:clear; clc;ezplot( '(x-2)A2+(y-3+2*x)A2-5' hold onezplot( '2*(x-3)A2+(y/3)A2-4' grid on;hold offf1=sym( '(x-2)A2+(y-3+2*x)A2=5' f2=sym( '2*(x-3)A2+(y/3)A2=4' x,y=solve(f1,f2, 'x' ,'y'); middle=x,y;intersection_x_y=double(middle)，-1,5,-8,8);，-1,5,-8,8)