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

中南大學(xué)MATLAB程序設(shè)計實(shí)踐材料科學(xué)與工程學(xué)院2023年3月26日一、編程實(shí)現(xiàn)“四階龍格－庫塔〔R-K〕方法求常微分方程〞，并舉一例應(yīng)用之?！緦?shí)例】采用龍格-庫塔法求微分方程：1、算法說明：在龍格-庫塔法中，四階龍格-庫塔法的局部截斷誤差約為o(h5)，被廣泛應(yīng)用于解微分方程的初值問題。其算法公式為：其中：2、流程圖：2.1、四階龍格－庫塔〔R-K〕方法流程圖：輸入待求微分方程、求解的自變量范圍、初值以及求解范圍內(nèi)的取點(diǎn)數(shù)等。輸入待求微分方程、求解的自變量范圍、初值以及求解范圍內(nèi)的取點(diǎn)數(shù)等。確定求解范圍內(nèi)的步長k=取點(diǎn)數(shù)？否求解：求解并輸出：是結(jié)束算法2.2、實(shí)例求解流程圖：輸入求解的自變量范圍輸入求解的自變量范圍求出待求簡單微分方程的真值解用MATLAB自帶函數(shù)ode23求解待求微分方程結(jié)束用自編函數(shù)四階龍格－庫塔〔R-K〕方法求解待求微分方程開始3、源程序代碼3.1、四階龍格－庫塔〔R-K〕方法源程序：function[x,y]=MyRunge_Kutta(fun,x0,xt,y0,PointNum,varargin)%Runge-Kutta方法解微分方程形為y'(t)=f(x,y(x))%此程序可解高階的微分方程。只要將其形式寫為上述微分方程的向量形式%函數(shù)f(x,y):fun%自變量的初值和終值：x0,xt%y0表示函數(shù)在x0處的值，輸入初值為列向量形式%自變量在[x0,xt]上取的點(diǎn)數(shù)：PointNum%varargin為可輸入項(xiàng)，可傳適當(dāng)參數(shù)給函數(shù)f(x,y)%x：所取的點(diǎn)的x值%y：對應(yīng)點(diǎn)上的函數(shù)值ifnargin<4|PointNum<=0PointNum=100;endifnargin<3y0=0;endy(1,:)=y0(:)';%初值存為行向量形式h=(xt-x0)/(PointNum-1);%計算步長x=x0+[0:(PointNum-1)]'*h;%得x向量值fork=1:(PointNum)%迭代計算f1=h*feval(fun,x(k),y(k,:),varargin{:});f1=f1(:)';%得公式k1f2=h*feval(fun,x(k)+h/2,y(k,:)+f1/2,varargin{:});f2=f2(:)';%得公式k2f3=h*feval(fun,x(k)+h/2,y(k,:)+f2/2,varargin{:});f3=f3(:)';%得公式k3f4=h*feval(fun,x(k)+h,y(k,:)+f3,varargin{:});f4=f4(:)';%得公式k4y(k+1,:)=y(k,:)+(f1+2*(f2+f3)+f4)/6;%得y(n+1)end3.2、實(shí)例求解源程序：%運(yùn)行四階R-K法clear,clc%去除內(nèi)存中的變量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構(gòu)造函數(shù)f(x,y)y0=0;%設(shè)定函數(shù)初值PointNum=5;%設(shè)定取點(diǎn)數(shù)[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)holdonplot(x1,y1,'bo',xr,yr,'r*')4、程序運(yùn)行結(jié)果：MyRunge_Kutta_x=00.50001.00001.50002.0000MyRunge_Kutta_y=00.39320.63180.77660.8645二、變成解決以下科學(xué)計算問題：〔一〕[例7-2-4]材料力學(xué)復(fù)雜應(yīng)力狀態(tài)的分析——Moore圓。1、程序說明：利用平面應(yīng)力狀態(tài)下斜截面應(yīng)力的一般公式，畫出任意平面應(yīng)力狀態(tài)下的應(yīng)力圓〔Moore圓〕，求出相應(yīng)平面應(yīng)力狀態(tài)下的主應(yīng)力〔、〕，并求出該應(yīng)力狀態(tài)下任意方位角的斜截面上的應(yīng)力、。2、程序流程圖：開始開始輸入待求應(yīng)力狀態(tài)的參數(shù)畫出應(yīng)力圓求某一方向角截面上的應(yīng)力?輸入方向角求出相應(yīng)正應(yīng)力、切應(yīng)力是否得出該應(yīng)力狀態(tài)下的主應(yīng)力求出主應(yīng)力平面方向角結(jié)束3、程序代碼：clear;clc;Sx=input('Sigma_x(MPa)=');%輸入該應(yīng)力狀態(tài)下的各應(yīng)力值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);%應(yīng)力圓一般方程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]);axisequal;%控制各坐標(biāo)軸的分度使其相等——使應(yīng)力圓變圓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);Tmin=min(Tau);b=axis;%提取坐標(biāo)軸邊界ps_array.Color='k';%控制坐標(biāo)軸顏色為黑色line([b(1),b(2)],[0,0],ps_array);%調(diào)整坐標(biāo)軸line([0,0],[b(3),b(4)],ps_array);b=axis;%提取坐標(biāo)軸邊界line([b(1),b(2)],[0,0],ps_array);%畫出x坐標(biāo)軸holdonplot(Sa,0,'.r')%標(biāo)出圓心text(Sa,0,'O');plot(Smax,0,'.r',Smin,0,'.r',Sa,Tmax,'.r',Sa,Tmin,'.r')%標(biāo)出最大、最小拉應(yīng)力、切應(yīng)力點(diǎn)text(Smax,0,'C');text(Smin,0,'D');text(Sa,Tmax,'E');text(Sa,Tmin,'F');%-----------此局部求某一斜截面上的應(yīng)力---------------------------------------------t=input('假設(shè)需求某一截面上的應(yīng)力，請輸入1；假設(shè)不求應(yīng)力，請輸入0：');whilet~=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('假設(shè)還需求其他截面上的應(yīng)力，請輸入1；假設(shè)要退出，請輸入0：');endholdoff%----------此局部求該應(yīng)力狀態(tài)下的主應(yīng)力--------------------------------------------Sigma_Max=SmaxSigma_Min=SminTau_Max=TmaxTau_Min=TminSigma1=Smax%得出主應(yīng)力Sigma3=Sminl=Sx-Sa;h=Txy;ratio=abs(h/l);%求主應(yīng)力平面方向角'主應(yīng)力平面方向角〔角度〕：'alpha_0=atan(ratio)/2*180/pi4、程序運(yùn)行結(jié)果：〔以為例〕Sigma_x(MPa)=100Sigma_y(MPa)=30Tau_xy(MPa)=-20假設(shè)需求某一截面上的應(yīng)力，請輸入1；假設(shè)不求應(yīng)力，請輸入0：1給出斜截面方向角：alpha=〔角度〕：30sigma=99.8205tau=20.3109假設(shè)還需求其他截面上的應(yīng)力，請輸入1；假設(shè)要退出，請輸入0：0Sigma_Max=105.3087Sigma_Min=24.6970Tau_Max=40.3109Tau_Min=-40.2963Sigma1=105.3087Sigma3=24.6970ans=主應(yīng)力平面方向角〔角度〕：alpha_0=14.8724〔二〕實(shí)驗(yàn)5〔橢圓的交點(diǎn)〕兩個橢圓可能具有0~4個交點(diǎn)，求以下兩個橢圓的所有交點(diǎn)坐標(biāo)：(1);(2)1、算法說明：此題相當(dāng)于求兩一個二元二次方程組的解，故為便于清楚地顯示出兩橢圓的相對位置，用ezplot函數(shù)把兩個橢圓畫在同一個坐標(biāo)系中，然后利用fsolve函數(shù)解方程組得到兩橢圓的交點(diǎn)即方程組的解。2、程序流程圖：開始開始依照所給方程畫出兩個橢圓用符號解法求兩橢圓的交點(diǎn)將符號解轉(zhuǎn)換成數(shù)值解結(jié)束3、程序代碼：clear;clc;ezplot('(x-2)^2+(y-3+2*x)^2-5',[-1,5,-8,8]);%畫第一個橢圓holdonezplot('2*(x-3)^2+(y/3)^2-4',[-1,5,-8,8]);%畫第二個橢圓gridon;%顯示網(wǎng)格holdofff1=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');%聯(lián)立兩個橢圓方程求解交點(diǎn)middle=[x,y];