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

1、欧拉近似方法求常微分方程朱翼1、编程实现以下科学计算算法，并举一例应用之。 “欧拉近似方法求常微分方程”算法说明：欧拉法是简单有效的常微分方程数值解法，欧拉法有多种形式的算法，其中简单欧拉法是一种单步递推算法。其基本原理为对简单的一阶方程的初值问题： y=f(x,y)其中 y(x0 )=y0欧拉法等同于将函数微分转换为数值微分，由欧拉公式可得 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 n

2、argin<6,trace=0;endt=t0;hmax=(tfinal-t)/16;1 / 19h=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) %主循环 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.' end if trace home,t,h,y end % Update the step si

4、ze if delta=0.0 h=min(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;tou

5、t=zeros(chunk,1);yout=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

6、=t+h; y=y+h*f; k=k+l;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。解题思路：首先建立例子所给函数的函数文件，调用函数my

7、euler，利用程序求解方程。将欧拉法解和精确解比较，由方程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('欧拉法','精

8、确解')例题流程图：f=f(x,y) f=-y+x+1;开始调用函数myeuler ， x,y=myeuler(f,0,1,1);求出结果y1=x+exp(-x)调用函数plot 绘图，比较 结束2、编程解决以下科学计算问题（1） 解题思路：建模： 材料力学中从弯矩求转角要经过一次不定积分, 而从转角求挠度又要经过一次不定积分, 在MATLAB中这却是非常简单的问题.可用cumsum函数作近似的不定积分,还可用更精确的函数cumtrapz来做不定积分。本题用cumsum函数来做.解题的关键还是在于正确地列写弯矩方程。本例中弯矩为程序代码：>> L=1; P=1000; L1

9、=1;%给出常数E = 200*109; I=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.

10、0098 -0.0122 -0.0146 -0.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.06

11、93 -0.0710 -0.0726 -0.0742 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

12、 -0.1087 Columns 64 through 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 9

13、1 through 99 -0.1251 -0.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

14、 -0.0323 -0.0362 -0.0403 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 -

15、0.2208 Columns 46 through 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

16、through 81 -0.5132 -0.5249 -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.84

17、59 -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 % 绘弯矩图流程图开始 L=1; P=1000; L1=1;E=200*109;I=2*10-5;调用linspace函数，将L分成100段n1=L1/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)*d

18、x结束（2）计算积分 解题思路：exp（-x2）是不可积函数，matlab中int积分显示不出积分结果，用到vpa函数控制其结果精度，表示出来。程序：>> syms x>> t=vpa(int (exp(-x.2)./(1+x.2),-inf,+inf),5)结果：t =1.3433解题思路：先建立内置函数，然后调用quad函数求积分。程序：>> fun=(x)tan(x)./x.(0.7);>> quad('fun',0,1)结果：ans = 0.9063解题思路：先建立内置函数，然后调用quad函数求积分。程序：>> fun=inline('exp(x)./(1-x.2).0.5');>> quad('fun',0,1)结果：ans =3.1044解题思路：这是一个二重积分，一般的方法是把二重积分化为二次积分，但由于二次积分命令int(int (1+x+y