龍格-庫(kù)塔法求微分方程2

1、MATLAB程序設(shè)計(jì)實(shí)踐課程考核一、編程實(shí)現(xiàn)“四階龍格庫(kù)塔（R-K）方法求常微分方程”，并舉一例應(yīng)用之?！緦?shí)例】采用龍格-庫(kù)塔法求微分方程：1、算法說(shuō)明：在龍格-庫(kù)塔法中，四階龍格-庫(kù)塔法的局部截?cái)嗾`差約為o(h5)，被廣泛應(yīng)用于解微分方程的初值問(wèn)題。其算法公式為：其中：2、流程圖：2.1、四階龍格庫(kù)塔（R-K）方法流程圖：輸入待求微分方程、求解的自變量范圍、初值以及求解范圍內(nèi)的取點(diǎn)數(shù)等。確定求解范圍內(nèi)的步長(zhǎng)k = 取點(diǎn)數(shù)？否求解：求解并輸出：是結(jié)束算法2.2、實(shí)例求解流程圖：輸入求解的自變量范圍求出待求簡(jiǎn)單微分方程的真值解用MATLAB自帶函數(shù)ode23求解待求微分方程結(jié)束用自編函數(shù)四階龍格

2、庫(kù)塔（R-K）方法求解待求微分方程開(kāi)始3、源程序代碼3.1、四階龍格庫(kù)塔（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

3、值%y：對(duì)應(yīng)點(diǎn)上的函數(shù)值 if nargin<4 | PointNum<=0 PointNum=100;endif nargin<3 y0=0;endy(1,:)=y0(:)' %初值存為行向量形式h=(xt-x0)/(PointNum-1); %計(jì)算步長(zhǎng) x=x0+0:(PointNum-1)'*h; %得x向量值for k=1:(PointNum) %迭代計(jì)算 f1=h*feval(fun,x(k),y(k,:),varargin:); f1=f1(:)' %得公式k1 f2=h*feval(fun,x(k)+h/2,y(k,:)+f1/2,var

4、argin:); 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)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*

5、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-')l

6、egend('真值','ode23','Rung-Kutta法解',2)hold onplot(x1,y1,'bo',xr,yr,'r*')4、程序運(yùn)行結(jié)果：MyRunge_Kutta_x = 0 0.5000 1.0000 1.5000 2.0000MyRunge_Kutta_y = 0 0.3932 0.6318 0.7766 0.8645二、編程解決以下科學(xué)計(jì)算問(wèn)題：（一）例7-2-4 材料力學(xué)復(fù)雜應(yīng)力狀態(tài)的分析Moore圓。1、程序說(shuō)明：利用平面應(yīng)力狀態(tài)下斜截面應(yīng)力的一般公式，畫出任意平面應(yīng)力狀態(tài)下的應(yīng)力圓

7、（Moore圓），求出相應(yīng)平面應(yīng)力狀態(tài)下的主應(yīng)力（、），并求出該應(yīng)力狀態(tài)下任意方位角的斜截面上的應(yīng)力、。2、程序流程圖：開(kāi)始輸入待求應(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=(

8、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);axis equal; %控制各坐標(biāo)軸的分度使其相等使應(yīng)力圓變圓title('應(yīng)力圓');xlabel('正應(yīng)力（MPa）');ylabel('剪應(yīng)力（MPa）');

9、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)軸hold onplot(Sa,0

10、,'.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('若需求某一截面上的應(yīng)力，請(qǐng)輸入1；若不求應(yīng)力

11、，請(qǐng)輸入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('若還需求其他截面上的應(yīng)力，請(qǐng)輸入1；若要退出，請(qǐng)輸入0：');endhold off%-此部分求該應(yīng)力狀態(tài)下的主應(yīng)力-Sigma_Max=SmaxSi

12、gma_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若需求某一截面上的應(yīng)力，請(qǐng)輸入1；若不求應(yīng)力，請(qǐng)輸入0：1給出斜截面方向角：alpha=（角度）：30sigma = 99.8205tau = 20.3109若還需求其

13、他截面上的應(yīng)力，請(qǐng)輸入1；若要退出，請(qǐng)輸入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)） 兩個(gè)橢圓可能具有0 4個(gè)交點(diǎn)，求下列兩個(gè)橢圓的所有交點(diǎn)坐標(biāo)：(1) ; (2) 1、算法說(shuō)明：此題相當(dāng)于求兩個(gè)二元二次方程組的解，故為便于清楚地顯示出兩橢圓的相對(duì)位置，用ezplot函數(shù)把兩個(gè)橢圓畫在同一個(gè)坐標(biāo)系中，然后利用fsolve函數(shù)解方程組得到兩橢圓的交點(diǎn)即方程組的解。2、程序流程圖：開(kāi)始依照所給方程畫出兩個(gè)橢圓用符號(hào)解法求兩橢圓的交點(diǎn)將符號(hào)解轉(zhuǎn)換成數(shù)值解結(jié)束3、程序代碼：clear; clc;ezplot('(x-2)2+(y-3+2*x)2-5',-1,5, -8,8); %畫第一個(gè)橢圓hold onezplot('2*(x-3)2+(y/3)2-4',-1,5, -8,8); %畫第二個(gè)橢圓grid on; %顯示網(wǎng)格hold