1、歐拉近似方法求常微分方程朱翼1、編程實現以下科學計算算法，并舉一例應用之?！皻W拉近似方法求常微分方程”算法說明：歐拉法是簡單有效的常微分方程數值解法，歐拉法有多種形式的算法，其中簡單歐拉法是一種單步遞推算法。其基本原理為對簡單的一階方程的初值問題：y'=f(x,y)其中y(x0)=y0歐拉法等同于將函數微分轉換為數值微分，由歐拉公式可得yn+1=yn+hf(xn,yn)程序代碼：functiontout,yout=myeuler(ypfun,t0,tfinal,y0,tol,trace)制始化pow=1/3;ifnargin<5,tol=1.e-3;endifnargin<

2、6,trace=0;endt=t0;hmax=(tfinal-t)/16;h=hmax/8;y=y0(:);chunk=128;tout=zeros(chunk,1);yout=zeros(chunk,length(y);k=1;tout(k)=t;yout(k,:)=y.'iftrace底圖clc,t,h,yendwhile(t<tfinal)&(t+h>t)%£循環ift+h>tfinal,h=tfinal-t;end%Computetheslopesf=feval(ypfun,t,y);f=f(:);%古計誤差并設定可接受誤差delta=nor

3、m(h*f,'inf');tau=tol*max(norm(y,'inf),1.0);哨誤差可接受時重寫解ifdelta<=taut=t+h;y=y+h*f;k=k+1;ifk>length(tout)tout=tout;zeros(chunk,1);yout=yout;zeros(chunk,length(y);endtout(k尸t;yout(k,:)=y.'endiftracehome,t,h,yend%Updatethestepsizeifdelta=0.0h=min(hmax,0.9*h*(tau/delta)Apow);endendif(

4、t<tfinal)dish('Singularitylikely.')tendtout=tout(1:k);yout=yout(1:k,:);流程圖:t+h>tfinalh=tfinal-t舉例：用歐拉法求y=-y+x+1,y(0)=1。解題思路：首先建立例子所給函數的函數文件，調用函數myeuler,利用程序求解方程。將歐拉法解和精確解比較，由方程f=-y+x+1可得到其精確解為y(x)=x+exp(-x)最后在同一圖窗中分別畫出兩圖。程序代碼：f.mfunctionf=f(x,y)f=-y+x+1;> >x,y=myeuler('f,0,1,

5、1);%J用程序求解方程> >y1=x+exp(-x);%方程f=-y+x+1的精確解> >plot(x,y,'-b',x,y1,'-r')%在同一圖窗將歐拉法解和精確解的圖畫出> >legend(歐拉法','精確解')例題流程圖:2、編程解決以下科學計算問題(1)L例72-£長為L的最胃梁加田723,在離固定魏L處事加力P.求它的轉角和提鹿.段梁矩為dx.EJE=200產UN/十和J2xTm,為巳旬.解題思路：建模：材料力學中從彎矩求轉角要經過一次不定積分，而從轉角求撓度又要經過一次不定積分

6、，在MATLAB中這卻是非常簡單的問題.可用cumsum函數作近似的不定積分，還可用更精確的函數cumtrapz來做不定積分。本題用cumsum函數來做.解題的關鍵還是在于正確地列寫彎矩方程。本例中彎(o<x<4)(Z:<x<Z)撓度Y=4dx程序代碼:>>L=1;P=1000;L1=1;%出常數E=200*10八9;1=2*10八-5;x=linspace(0,L,101);dx=L/100;%Fx分100段n1=L7dx+1;%確定x=L1處對應的下標M1=-P*(L1-x(1:n1);%第一段彎矩賦值M2=zeros(1,101-n1);%第二段彎矩賦

7、值(全為零)M=M1,M2;%全梁的彎矩A=cumsum(M)*dx/(E*I)%對彎矩積分求轉角Y=cumsum(A)*dx%對轉角積分求撓度A=1.0e-003*-0.0146Columns1through9-0.0025-0.0050-0.0074-0.0098-0.0122-0.0170-0.0193-0.0216Columns10through18-0.0239-0.0261-0.0283-0.0305-0.0327-0.0349-0.0370-0.0391-0.0412Columns19through27-0.0432-0.0452-0.0472-0.0492-0.0512-0.0

8、531-0.0550-0.0569-0.0587Columns28through36-0.0605-0.0623-0.0641-0.0659-0.0676-0.0693-0.0710-0.0726-0.0742Columns37through45-0.0759-0.0774-0.0790-0.0805-0.0820-0.0835-0.0849-0.0863-0.0877-0.0944-0.0956Columns46through54-0.0891-0.0905-0.0918-0.0931-0.0968-0.0980-0.0992Columns55through63-0.1004-0.1015-

9、0.1026-0.1037-0.1047-0.1057-0.1067-0.1077-0.1087Columns64through72-0.1096-0.1105-0.1114-0.1122-0.1130-0.1138-0.1146-0.1154-0.1161Columns73through81-0.1168-0.1175-0.1181-0.1187-0.1194-0.1199-0.1205-0.1210-0.1215Columns82through90-0.1220-0.1224-0.1229-0.1232-0.1236-0.1240-0.1243-0.1246-0.1249Columns91

10、through99-0.1251-0.1253-0.1255-0.1257-0.1259-0.1260-0.1261-0.1262-0.1262Columns100through101-0.1262-0.12621.0e-004*-0.0025-0.0037-0.0052-0.0218-0.0251-0.0286Columns1through9-0.0002-0.0007-0.00150.0069-0.0088-0.0109Columns10through18-0.0133-0.0159-0.01880.0323-0.0362-0.0403Columns19through27-0.0446-0

11、.0492-0.0539-0.0588-0.0639-0.0692-0.0747-0.0804-0.0863Columns28through36-0.0924-0.0986-0.1050-0.1116-0.1184-0.1253-0.1324-0.1397-0.1471Columns37through45-0.1547-0.1624-0.1703-0.1783-0.1865-0.1949-0.2034-0.2120-0.2208Columns46through54-0.2297-0.2388-0.2479-0.2572-0.2667-0.2762-0.2859-0.2957-0.3057Col

12、umns55through63- 0.3157-0.3258-0.3361-0.3782-0.3890-0.3998Columns64through72- 0.4108-0.4218-0.4330-0.4784-0.4899-0.5015Columns73through81- 0.5132-0.5249-0.5367-0.5846-0.5967-0.6088Columns82through90- 0.6210-0.6333-0.6456-0.6951-0.7075-0.7200Columns91through99- 0.7325-0.7451-0.7576-0.3465-0.3569-0.36

13、75-0.4442-0.4555-0.4669-0.5486-0.5606-0.5726-0.6579-0.6703-0.6827-0.7702-0.7828-0.7954-0.8080-0.8206-0.8332-0.8459-0.8585>>subplot(3,1,1),plot(x,M),grid%繪彎矩圖subplot(3,1,2),plot(x,A),grid%繪彎矩圖subplot(3,1,3),plot(x,Y),grid%繪彎矩圖EditTextGoCellTo.olsDebugDezktopWindowHelpp31同£*：+ffleg-1.0+.11x

14、喊碰GtL-L=l;P=1000;Ll=l;蜓合出常數2-E=200*10*9；1=2*10*-5；3 -x=lizisp&ce(0L,101);dx=L/10口；鬻乂分1口口段4 -nl=Ll/dM+l;%確定k=LI,界對向的下標5 -N1=_F*(Ll-s(1;%第一段彎矩賦值6 -M2=zeros(L101-nlk%第二段彎矩呢:值(全為雪)7 -M=MLM2;抬全梁的彎矩S-A=curnwumCE*工)%對彎矩枳分求轉角勺-Y-匚umwum(A)*dx%對轉角枳分求撓度10 -subplotC3flDplot(je,grid找繪彎矩圖11 -與ubpl”1,2)11PLota

15、,A)/grid%筵彎矩圖12 -subplot0,L3),Plot(x,Y),grid抬繪彎矩圖流程圖n1=L"dx+1;M1=-P*(L1-x(1:n1);M2=zeros(1,101-n1);M=M1,M2;調用cumsum函數求A=cumsum(M)*dx/(E*I)Y=cumsum(A)*dx結束(2)計算積分1一“解題思路：exp(-xA2)是不可積函數，matlab中int積分顯示不出積分結果，用到vpa函數控制其結果精度，表示出來。程序：>>symsx>>t=vpa(int(exp(-x.A2)./(1+x.A2),-inf,+inf),5)結

16、果：t=1.3433J1工解題思路：先建立內置函數，然后調用quad函數求積分。程序：>>fun=(x)tan(x)./x.A(0.7);>>quad('fun',0,1)結果：ans=0.9063(普，L;解題思路：先建立內置函數，然后調用quad函數求積分。程序：>>fun=inline('exp(x)./(1-x.A2).A0.5');>>quad('fun',0,1)結果：ans=3.1044“(l+彳+/6dm為彳口+/<2x解題思路：這是一個二重積分，一般的方法是把二重積分化為二次積分，但由于二次積分命