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

1、欧拉近似方法求常微分方程朱翼1、编程实现以下科学计算算法，并举一例应用之。“欧拉近似方法求常微分方程”算法说明：欧拉法是简单有效的常微分方程数值解法，欧拉法有多种形式的算法，其中简单欧拉法是一种单步递推算法。其基本原理为对简单的一阶方程的初值问题：y-f(x,y)其中 y(x0 )=y()欧拉法等同于将函数微分转换为数值微分，由欧拉公式可得yn+1 =y n+hf(x n ,y n)程序代码：function tout,yout=myeuler(ypfun,t0,tfinal,y0,tol,trace)%J始化pow=1/3;if nargin<5,tol=1.e-3; endif na

2、rgin<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)%£循环if t+h>tfinal,h=tfinal-t; end% Compute the slopesf=feval(ypfun,t,y);f=f(:);%古计误差并

3、设定可接受误差delta=norm(h*f, 'inf');tau=tol*max(norm(y,'inf),1.0);%i误差可接受时重写解if delta<=taut=t+h;y=y+h*f;k=k+1;if k>length(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.9*h

4、*(tau/delta)Apow);endendif (t<tfinal)dish( 'Singularity likely.' )tendtout=tout(1:k);yout=yout(1:k,:);流程图:h=tfinal-th=min(hmax,0.9*h*(tau/delta)Apow)举例：用欧拉法求 y=-y+x+1 , y (0) =1。解题思路：首先建立例子所给函数的函数文件，调用函数 myeuler,利用程序 求解方程。将欧拉法解和精确解比较，由方程f=-y+x+1可得到其精确 解为y (x) =x+exp(-x)最后在同一图窗中分别画出两图。程序代码

5、：f. mfunction f=f(x,y)f=-y+x+1;>>x,y=myeuler(T,0,1,1);%用程序求解方程>> y1=x+exp(-x);%方程 f=-y+x+1 的精确解>>plot(x,y,'-b',x,y1,'-r') %在同一图窗将欧拉法解和精确解的图画出>> legend(欧拉法','精确解')例题流程图:2、编程解决以下科学计算问题(1)长为L的揪曾梁如为723，左端团,在离固定端L业澹加力P.术它的的角加援度.设梁E = 20Q，70'N/m'

6、和为巳知,解题思路：建模：材料力学中从弯矩求转角要经过一次不定积分, 而从转角求挠度又要经过一次不定积分，在MATLAB中这却 是非常简单的问题.可用cumsum函数作近似的不定积分，还可 用更精确的函数cumtrapz来做不定积分。本题用cumsum函 数来做.解题的关键还是在于正确地列写弯矩方程。本例中弯矩为转角 A = dx,Jo EJ(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/

7、100;%F x 分 100 段n1=L7dx+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.0170-0.0193-0.0216Columns 10 through 18-0.023

8、9-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.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.0835-0.0849-0.08

9、63-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.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.1

10、181-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-0.1246-0.1249Columns 91 through 99-0.1251-0.1253-0.1255-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.0069-0.0088-0.0

11、109Columns 10 through 18-0.0133-0.0159-0.0188-0.1257-0.1259-0.1260-0.0025-0.0037-0.0052-0.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.1471Col

12、umns 37 through 45-0.1547-0.1624-0.1703-0.1783-0.1865-0.1949-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.3998-0.4555-0.4669Columns 64 through 72-0.4108-0.4218-0

13、.4330-0.4442-0.4784-0.4899-0.5015-0.5726Columns 73 through 81-0.5132-0.5249-0.5367-0.5486-0.5606-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 thr

14、ough 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% 绘弯矩图Fil Edi t Text Go Cell To.ols Debug Dezktop Window Help：Cl 诰唱 ©百昌系冒匚冒 "d+ |_ x 蜷癖GtL - L=lf P=1000; Ll = l ：艇出常数2 -E = 200*10*9； 1=2*10"-5；3 -x = li

15、zisp&ce (0 L, 101 ) ; dxL/lOO;段4 - nl=Ll/dM+l;%晞定试二L1炫对成的下标5 - Nl = -F* ( Ll-s (1;% 第一段穹矩嗽值6 - M2 = zeros (L 101-nlk%第二段弯矩赋值(全为耍)7 - M = ML M2;肾全梁的弯矩S -A = cumsum(M)*dx/ CE*I) %对弯矩积金'求转角勺-Y - cumwum (A)*dx%对转角积分求撬度10 -subplot C3f lDplot(je,gri d帘绘弯矩图11 -subplot (3 I, 2)/ plot6c,A)gri d% 绘弯矩图

16、12 -subplot 0, L 3), Plot(x,Y), gri d需绘弯矩匿®8n1=LTdx+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)计算积分 J - 解题思路：exp (-xA2)是不可积函数，matlab中int积分显示不出积分结果，用到vpa函数控制其结果精度，表示出来。程序：>> syms x >> t=vpa(int (exp(-x.A2)./(1+x.A2),-inf,

17、+inf),5)结果：t =1.3433-1 T解题思路：先建立内置函数，然后调用 quad函数求积分。程序：>> fun=(x)tan(x)./x.A(0.7);>> quad('fun',0,1)结果：ans = 0.9063，",解题思路：先建立内置函数，然后调用 quad函数求积分。程序：>> fun=inline('exp(x)./(1-x.A2).A0.5');>> quad('fun',0,1)结果：ans =3.1044|(1+了 + 物以£)为了+尸<2x解题思路：这是一个二重积分，一般的方法是把二重积分化为二次积