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

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

2、gin<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.'if trace %繪圖 clc,t,h,yendwhile (t<tfinal)&(t+h>t) %主循環(huán)if t+h>tfinal,h=tfinal-t;end% Compute the slopes f=feval(ypfun,t,y);f=f(:);%估計誤差并設定可

3、接受誤差 delta=norm(h*f,'inf'); tau=tol*max(norm(y,'inf'),1.0);%當誤差可接受時重寫解if delta<=tau t=t+h; y=y+h*f; k=k+1;if k>length(tout) tout=tout;zeros(chunk,1); yout=yout;zeros(chunk,length(y);end tout(k)=t; yout(k,:)=y.'endif trace home,t,h,yend% Update the step sizeif delta=0.0 h=mi

4、n(hmax,0.9*h*(tau/delta)pow);endendif (t<tfinal) dish('Singularity likely.') tendtout=tout(1:k);yout=yout(1:k,:);流程圖：開始tout,yout=myeuler(ypfun,t0,tfinal,y0,tol,trace)Pow=1/3Ynargin<5Ntol=1.e-3Ynargin<6trace=0t=t0;hmax=(tfinal-t)/16;h=hmax/8;y=y0(:);chunk=128;tout=zeros(chunk,1);yout

5、=zeros(chunk,length(y); k=1;tout(k)=t;yout(k,:)=y.'NYtrace= 1N輸出t,h,yNt<tfinal & t+h>tYYt+h>tfinalh=tfinal-tNf=feval(ypfun,t,y); f=f(:); delta=norm(h*f,'inf'); tau=tol*max(norm(y,'inf')1.0); tau=tol*max(norm(y,'inf')1.0);Ydelta <=tauNt=t+h; y=y+h*f; k=k+l;

6、k>length(tout)tout=tout;zeros(chunk,1); yout=yout;zeros(chunk,length(y);tout(k)=t; yout(k,:)=y,'Ytrace=1輸出home ,t ,h,yNYdelta=0.0 h=min(hmax,0.9*h*(tau/delta)pow)NYt<tfinal輸出 disp tN輸出tout=tout(l:k);yout=yout(l:k,:);結束舉例：用歐拉法求y=-y+x+1，y（0）=1。解題思路：首先建立例子所給函數(shù)的函數(shù)文件，調(diào)用函數(shù)myeuler，利用程序求解方程。將歐拉法解和

7、精確解比較，由方程f=-y+x+1可得到其精確解為y（x）=x+exp(-x)。最后在同一圖窗中分別畫出兩圖。程序代碼：fmfunction f=f(x,y)f=-y+x+1;>>x,y=myeuler('f',0,1,1); %利用程序求解方程>> y1=x+exp(-x); %方程f=-y+x+1的精確解>>plot(x,y,'-b',x,y1,'-r')%在同一圖窗將歐拉法解和精確解的圖畫出>> legend('歐拉法','精確解')例題流程圖：f=f(x,y)

8、 f=-y+x+1;開始調(diào)用函數(shù)myeuler ， x,y=myeuler(f,0,1,1);求出結果y1=x+exp(-x)調(diào)用函數(shù)plot 繪圖，比較 結束2、編程解決以下科學計算問題（1）解題思路：建模： 材料力學中從彎矩求轉角要經(jīng)過一次不定積分, 而從轉角求撓度又要經(jīng)過一次不定積分, 在MATLAB中這卻是非常簡單的問題.可用cumsum函數(shù)作近似的不定積分,還可用更精確的函數(shù)cumtrapz來做不定積分。本題用cumsum函數(shù)來做.解題的關鍵還是在于正確地列寫彎矩方程。本例中彎矩為程序代碼：>> L=1; P=1000; L1=1;%給出常數(shù)E = 200*109; I=

9、2*10-5;x = linspace(0,L,101); dx=L/100;%將x分100段n1=L1/dx+1;% 確定x=L1處對應的下標M1 = -P*( L1-x(1:n1); % 第一段彎矩賦值M2 = zeros(1,101-n1); % 第二段彎矩賦值（全為零）M = M1,M2;% 全梁的彎矩A = cumsum(M)*dx/(E*I)% 對彎矩積分求轉角Y = cumsum(A)*dx % 對轉角積分求撓度A = 1.0e-003 * Columns 1 through 9 -0.0025 -0.0050 -0.0074 -0.0098 -0.0122 -0.0146 -0

10、.0170 -0.0193 -0.0216 Columns 10 through 18 -0.0239 -0.0261 -0.0283 -0.0305 -0.0327 -0.0349 -0.0370 -0.0391 -0.0412 Columns 19 through 27 -0.0432 -0.0452 -0.0472 -0.0492 -0.0512 -0.0531 -0.0550 -0.0569 -0.0587 Columns 28 through 36 -0.0605 -0.0623 -0.0641 -0.0659 -0.0676 -0.0693 -0.0710 -0.0726 -0.0

11、742 Columns 37 through 45 -0.0759 -0.0774 -0.0790 -0.0805 -0.0820 -0.0835 -0.0849 -0.0863 -0.0877 Columns 46 through 54 -0.0891 -0.0905 -0.0918 -0.0931 -0.0944 -0.0956 -0.0968 -0.0980 -0.0992 Columns 55 through 63 -0.1004 -0.1015 -0.1026 -0.1037 -0.1047 -0.1057 -0.1067 -0.1077 -0.1087 Columns 64 thr

12、ough 72 -0.1096 -0.1105 -0.1114 -0.1122 -0.1130 -0.1138 -0.1146 -0.1154 -0.1161 Columns 73 through 81 -0.1168 -0.1175 -0.1181 -0.1187 -0.1194 -0.1199 -0.1205 -0.1210 -0.1215 Columns 82 through 90 -0.1220 -0.1224 -0.1229 -0.1232 -0.1236 -0.1240 -0.1243 -0.1246 -0.1249 Columns 91 through 99 -0.1251 -0

13、.1253 -0.1255 -0.1257 -0.1259 -0.1260 -0.1261 -0.1262 -0.1262 Columns 100 through 101 -0.1262 -0.1262Y = 1.0e-004 * Columns 1 through 9 -0.0002 -0.0007 -0.0015 -0.0025 -0.0037 -0.0052 -0.0069 -0.0088 -0.0109 Columns 10 through 18 -0.0133 -0.0159 -0.0188 -0.0218 -0.0251 -0.0286 -0.0323 -0.0362 -0.040

14、3 Columns 19 through 27 -0.0446 -0.0492 -0.0539 -0.0588 -0.0639 -0.0692 -0.0747 -0.0804 -0.0863 Columns 28 through 36 -0.0924 -0.0986 -0.1050 -0.1116 -0.1184 -0.1253 -0.1324 -0.1397 -0.1471 Columns 37 through 45 -0.1547 -0.1624 -0.1703 -0.1783 -0.1865 -0.1949 -0.2034 -0.2120 -0.2208 Columns 46 throu

15、gh 54 -0.2297 -0.2388 -0.2479 -0.2572 -0.2667 -0.2762 -0.2859 -0.2957 -0.3057 Columns 55 through 63 -0.3157 -0.3258 -0.3361 -0.3465 -0.3569 -0.3675 -0.3782 -0.3890 -0.3998 Columns 64 through 72 -0.4108 -0.4218 -0.4330 -0.4442 -0.4555 -0.4669 -0.4784 -0.4899 -0.5015 Columns 73 through 81 -0.5132 -0.5

16、249 -0.5367 -0.5486 -0.5606 -0.5726 -0.5846 -0.5967 -0.6088 Columns 82 through 90 -0.6210 -0.6333 -0.6456 -0.6579 -0.6703 -0.6827 -0.6951 -0.7075 -0.7200 Columns 91 through 99 -0.7325 -0.7451 -0.7576 -0.7702 -0.7828 -0.7954 -0.8080 -0.8206 -0.8332 Columns 100 through 101 -0.8459 -0.8585>> subp

17、lot(3,1,1),plot(x,M),grid % 繪彎矩圖subplot(3,1,2),plot(x,A),grid % 繪彎矩圖subplot(3,1,3),plot(x,Y),grid % 繪彎矩圖流程圖開始L=1; P=1000; L1=1;E=200*109;I=2*10-5;調(diào)用linspace函數(shù)，將L分成100段n1=L1/dx+1; M1 = -P*( L1-x(1:n1); M2 = zeros(1,101-n1);M= M1,M2;調(diào)用cumsum函數(shù)求A=cumsum(M)*dx/(E*I) Y=cumsum(A)*dx結束（2）計算積分解題思路：exp（-x2）是

18、不可積函數(shù)，matlab中int積分顯示不出積分結果，用到vpa函數(shù)控制其結果精度，表示出來。程序：>> syms x>> t=vpa(int (exp(-x.2)./(1+x.2),-inf,+inf),5)結果：t =1.3433解題思路：先建立內(nèi)置函數(shù)，然后調(diào)用quad函數(shù)求積分。程序：>> fun=(x)tan(x)./x.(0.7);>> quad('fun',0,1)結果：ans = 0.9063解題思路：先建立內(nèi)置函數(shù)，然后調(diào)用quad函數(shù)求積分。程序：>> fun=inline('exp(x)./(1-x.2).0.5');>> quad('fun',0,1)結果：ans =3.1044解題思路：