歐拉近似方法求常微分方程

1、歐拉近似方法求常微分方程朱翼1、編程實(shí)現(xiàn)以下科學(xué)計(jì)算算法，并舉一例使用之?！皻W拉近似方法求常微分方程”算法說明： 歐拉法是簡(jiǎn)單有效的常微分方程數(shù)值解法， 歐拉法有多種形式的 算法，其中簡(jiǎn)單歐拉法是一種單步遞推算法。其基本原理為 對(duì)簡(jiǎn)單的一階方程的初值問題：y=f(x,y)其中 y(x0 )=y0歐拉法等同于將函數(shù)微分轉(zhuǎn)換為數(shù)值微分，由歐拉公式可得yn+1 =y n+hf(x n ,y n)程序代碼：function tout,yout=myeuler(ypfun,t0,tfinal,y0,tol,trace)%初始化pow=1/3;if nargin5,tol=1.e-3; endif nar

2、gin6,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.;%繪圖%主循環(huán)if traceclc,t,h,yendwhile (tt)if t+htfinal,h=tfinal-t; end % Compute the slopes f=feval(ypfun,t,y);f=f(:);%估計(jì)誤差并設(shè)定可接受誤差 delta=norm(h*f, inf ); tau=tol

3、*max(norm(y, inf ),1.0);%當(dāng)誤差可接受時(shí)重寫解if deltalength(tout)tout=tout;zeros(chunk,1);yout=yout;zeros(chunk,length(y); endtout(k)=t; yout(k,:)=y.; endif trace home,t,h,y end% Update the step sizeif delta=0.0 h=min(hmax,0.9*h*(tau/delta)pow); endendif (ttfinal)dish( Singularity likely. )tend tout=tout(1:k)

4、; yout=yout(1:k,:);流程圖：開始Po3YNNtrace=NNttfinalNYNklengt(tout)f)1.0)t+ht Ynodrmel(tya, in=ft)a1u.tout,yout=myeuler(ypfun,t0,tfinal,y0,tol,trace)Y t=t0;hmax=(tfinal-t)/16;h=hmax/8;y=y0(:);chunk=128;nargin6tout=zeros(chunk,1);yout=zeros(chunk,length(y);k=1;tout(k)=t;yout(k,:)=y.;narginx,y=myeuler(f,0,1

5、,1); %利用程序求解方程 y1=x+exp(-x);%方程 f=-y+x+1 的精確解plot(x,y,-b,x,y1,-r) %在同一圖窗將歐拉法解和精確解的圖畫出 legend(歐 拉法 ,精確解 )例題流程圖： 開始f=f(x,y)f=-y+x+1;調(diào)用函數(shù)myeuler ，x,y=myeuler( f ,0,1,1);求出y結(jié)1果=x+exp(-x)2、1)解題思路 ：建模： 材料力學(xué)中從彎矩求轉(zhuǎn)角要經(jīng)過一次不定積分 , 而從轉(zhuǎn)角求撓度又要經(jīng)過一次不定積分 , 在 MATLAB中這卻 是非常簡(jiǎn)單的問題 .可用 cumsum函數(shù)作近似的不定積分 ,還可 用更精確的函數(shù) cumtrap

6、z 來做不定積分。本題用 cumsum 函 數(shù)來做.解題的關(guān)鍵還是在于正確地列寫彎矩方程。本例中彎 矩為程序代碼： L=1; P=1000; L1=1; 給%出常數(shù)E = 200*109; I=2*10-5;x = linspace(0,L,101); dx=L/100;%將 x分 100段 n1=L1/dx+1; % 確定 x=L1 處對(duì)應(yīng)的下標(biāo) M1 = -P*( L1-x(1:n1); % 第一段彎矩賦值M2 = zeros(1,101-n1); % 第二段彎矩賦值(全為零)M = M1,M2;% 全梁的彎矩A = cumsum(M)*dx/(E*I)% 對(duì)彎矩積分求轉(zhuǎn)角Y = cums

7、um(A)*dx% 對(duì)轉(zhuǎn)角積分求撓度1.0e-003 *Columns 1 through 9-0.0025 -0.0050-0.0074-0.0098-0.0122-0.0146-0.0170 -0.0193 -0.0216Columns 10 through 18-0.0239 -0.0261-0.0283-0.0305-0.0327-0.0349-0.0370 -0.0391 -0.0412Columns 19 through 27-0.0432 -0.0452 -0.0472-0.0492-0.0512-0.0531-0.0550 -0.0569 -0.0587Columns 28 t

8、hrough 36-0.0605 -0.0623 -0.0641-0.0659-0.0676-0.0693-0.0710 -0.0726 -0.0742Columns 37 through 45-0.0759 -0.0774 -0.0790-0.0805-0.0820-0.08350.0849 -0.0863 -0.0877Columns 46 through 54-0.0891 -0.0905 -0.0918-0.0931-0.0944-0.09560.0968 -0.0980 -0.0992Columns 55 through 63-0.1004 -0.1015 -0.1026-0.103

9、7-0.1047-0.10570.1067 -0.1077 -0.1087Columns 64 through 72-0.1096 -0.1105 -0.1114-0.1122-0.1130-0.11380.1146 -0.1154 -0.1161Columns 73 through 81-0.1168 -0.1175 -0.1181-0.1187-0.1194-0.11990.1205 -0.1210 -0.1215Columns 82 through 90-0.1220 -0.1224 -0.1229-0.1232-0.1236-0.12400.1243 -0.1246 -0.1249Co

10、lumns 91 through 99-0.1251 -0.1253 -0.1255-0.1257-0.1259-0.1260-0.1261 -0.1262 -0.1262Columns 100 through 101-0.1262 -0.12621.0e-004 *Columns 1 through 9-0.0002 -0.0007 -0.0015-0.0025-0.0037-0.00520.0069 -0.0088 -0.0109Columns 10 through 18-0.0133 -0.0159 -0.0188-0.0218-0.0251-0.02860.0323 -0.0362 -

11、0.0403Columns 19 through 27-0.0446 -0.0492 -0.0539-0.0588-0.0639-0.06920.0747 -0.0804 -0.0863Columns 28 through 36-0.0924 -0.0986 -0.1050-0.1116-0.1184-0.12530.1324 -0.1397 -0.1471Columns 37 through 45-0.1547 -0.1624 -0.1703-0.1783-0.1865-0.19490.2034 -0.2120 -0.2208Columns 46 through 54-0.2297 -0.2

12、388 -0.2479-0.2572-0.2667-0.27620.2859 -0.2957 -0.3057Columns 55 through 63-0.3157 -0.3258 -0.3361-0.3465-0.3569-0.3675-0.3782 -0.3890 -0.3998Columns 64 through 72-0.4108 -0.4218 -0.4330-0.4442-0.4555-0.4669-0.4784 -0.4899 -0.5015Columns 73 through 81-0.5132 -0.5249 -0.5367-0.5486-0.5606-0.5726-0.58

13、46 -0.5967 -0.6088Columns 82 through 90-0.6210 -0.6333 -0.6456-0.6579-0.6703-0.6827-0.6951 -0.7075 -0.7200Columns 91 through 99-0.7325 -0.7451 -0.7576-0.7702-0.7828-0.7954-0.8080 -0.8206 -0.8332Columns 100 through 101-0.8459 -0.8585 subplot(3,1,1),plot(x,M),grid % 繪彎矩圖 subplot(3,1,2),plot(x,A),grid%

14、 繪彎矩圖subplot(3,1,3),plot(x,Y),grid% 繪彎矩圖流程圖開始解題思路 ：exp( -x2)是不可積函數(shù)， matlab 中 int 積分顯示不出積 分結(jié)果，用到 vpa 函數(shù)控制其結(jié)果精度，表示出來。程序： syms x t=vpa(int (exp(-x.2)./(1+x.2),-inf,+inf),5)結(jié)果：t =1.3433解題思路 ：先建立內(nèi)置函數(shù)，然后調(diào)用 quad函數(shù)求積分。程序： fun=(x)tan(x)./x.(0.7); quad(fun,0,1)結(jié)果：ans = 0.9063解題思路 ：先建立內(nèi)置函數(shù)，然后調(diào)用 quad函數(shù)求積分。程序： fun=inline(exp(x)./(1-x.2).0.5); quad(fun,0,1)結(jié)果：ans =3.1044解題思路 ：這是一個(gè)二重積分， 一般的方法是把二重積分化為二次積 分，但由于二次積分命令int(int (1+x+y.2,y,-sq