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 nargint) %主循环if t+htfinal,h=tfin
2、al-t;end% Compute the slopesf=feval(ypfun,t,y);f=f(:);%估计误差并设定可接受误差delta=norm(h*f,inf);tau=tol*max(norm(y,inf),1.0);%当误差可接受时重写解if deltalength(tout)tout=tout;zeros(chunk,1);yout=yout;zeros(chunk,length(y);endtout(k)=t;yout(k,:)=y.;endif tracehome,t,h,yend% Update the step sizeif delta=0.0h=min(hmax,0
3、.9*h*(tau/delta)pow);endendif (ttNt+htfinalYYh=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);delta length(tout)tout=tout;zeros(chunk,1); yout=yout;zeros(chunk,length(y);tout(k)=t; yout(k,:)=y,;trace=1输出 home ,t ,h,yYNdelta=0.0 Yh=min
4、(hmax,0.9*h*(tau/delta)pow)N举例：用欧拉法求y=-y+x+1，y（0）=1。解题思路：首先建立例子所给函数的函数文件，调用函数myeuler，利用程序求解方程。将欧拉法解和精确解比较，由方程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
5、(欧拉法,精确解)t L=1; P=1000; L1=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 thr
6、ough 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 through 36 -0.0605 -0.062
7、3 -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.0835 -0.0849 -0.0863 -0.0877Columns 46 through 54 -0.0891 -0.0905 -0.0918 -0.0931 -0.0944 -0.0956 -0.0968 -0.0980 -0.0992Columns 55 through 63 -0.1004 -0.1015 -0.1026 -0.1037 -0.
8、1047 -0.1057 -0.1067 -0.1077 -0.1087Columns 64 through 72 -0.1096 -0.1105 -0.1114 -0.1122 -0.1130 -0.1138 -0.1146 -0.1154 -0.1161Columns 73 through 81 -0.1168 -0.1175 -0.1181 -0.1187 -0.1194 -0.1199 -0.1205 -0.1210 -0.1215Columns 82 through 90 -0.1220 -0.1224 -0.1229 -0.1232 -0.1236 -0.1240 -0.1243
9、-0.1246 -0.1249Columns 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.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.0109Columns 10 through 18 -0.0133 -0.0159 -0.0188 -0.
10、0218 -0.0251 -0.0286 -0.0323 -0.0362 -0.0403Columns 19 through 27 -0.0446 -0.0492 -0.0539 -0.0588 -0.0639 -0.0692 -0.0747 -0.0804 -0.0863Columns 28 through 36 -0.0924 -0.0986 -0.1050 -0.1116 -0.1184 -0.1253 -0.1324 -0.1397 -0.1471Columns 37 through 45 -0.1547 -0.1624 -0.1703 -0.1783 -0.1865 -0.1949
11、-0.2034 -0.2120 -0.2208Columns 46 through 54 -0.2297 -0.2388 -0.2479 -0.2572 -0.2667 -0.2762 -0.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.50
12、15Columns 73 through 81 -0.5132 -0.5249 -0.5367 -0.5486 -0.5606 -0.5726 -0.5846 -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
13、 101 -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 % 绘弯矩图流程图开始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=cu
14、msum(A)*dx（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.
15、5); quad(fun,0,1)结果：ans =3.1044解题思路：这是一个二重积分，一般的方法是把二重积分化为二次积分，但由于二次积分命令int(int (1+x+y.2,y,-sqrt(2*x-x.2),sqrt(2*x-x.2),x,0,1)无法求出结果，故将其转化成极坐标。积分函数为r.*(1+r.*cos(a)+r.2.*(sin(a).2)其中 a 的范围：-pi/2 到 pi/2；r 的范围：0 到 2.*cos(a)。程序： syms a r int(int(r+r.2.*cos(a)+r.3.*(sin(a).2,r,0,2.*cos(a),a,-pi/2,pi/2)结果：ans =(9*pi)/4
