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

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

2、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.'iftrace底图clc,t,h,yendwhile(t<tfinal)&(t+h>t)%£循环ift+h>tfinal,h=tfinal-t;end%Computetheslopesf=feval(ypfun,t,y);f=f(:);%古计误差并设定可接受误差delta=nor

3、m(h*f,'inf');tau=tol*max(norm(y,'inf),1.0);哨误差可接受时重写解ifdelta<=taut=t+h;y=y+h*f;k=k+1;ifk>length(tout)tout=tout;zeros(chunk,1);yout=yout;zeros(chunk,length(y);endtout(k尸t;yout(k,:)=y.'endiftracehome,t,h,yend%Updatethestepsizeifdelta=0.0h=min(hmax,0.9*h*(tau/delta)Apow);endendif(

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

5、1);%J用程序求解方程> >y1=x+exp(-x);%方程f=-y+x+1的精确解> >plot(x,y,'-b',x,y1,'-r')%在同一图窗将欧拉法解和精确解的图画出> >legend(欧拉法','精确解')例题流程图:2、编程解决以下科学计算问题(1)L例72-£长为L的最胃梁加田723,在离固定魏L处事加力P.求它的转角和提鹿.段梁矩为dx.EJE=200产UN/十和J2xTm,为巳旬.解题思路：建模：材料力学中从弯矩求转角要经过一次不定积分，而从转角求挠度又要经过一次不定积分

6、，在MATLAB中这却是非常简单的问题.可用cumsum函数作近似的不定积分，还可用更精确的函数cumtrapz来做不定积分。本题用cumsum函数来做.解题的关键还是在于正确地列写弯矩方程。本例中弯(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/100;%Fx分100段n1=L7dx+1;%确定x=L1处对应的下标M1=-P*(L1-x(1:n1);%第一段弯矩赋值M2=zeros(1,101-n1);%第二段弯矩赋

7、值(全为零)M=M1,M2;%全梁的弯矩A=cumsum(M)*dx/(E*I)%对弯矩积分求转角Y=cumsum(A)*dx%对转角积分求挠度A=1.0e-003*-0.0146Columns1through9-0.0025-0.0050-0.0074-0.0098-0.0122-0.0170-0.0193-0.0216Columns10through18-0.0239-0.0261-0.0283-0.0305-0.0327-0.0349-0.0370-0.0391-0.0412Columns19through27-0.0432-0.0452-0.0472-0.0492-0.0512-0.0

8、531-0.0550-0.0569-0.0587Columns28through36-0.0605-0.0623-0.0641-0.0659-0.0676-0.0693-0.0710-0.0726-0.0742Columns37through45-0.0759-0.0774-0.0790-0.0805-0.0820-0.0835-0.0849-0.0863-0.0877-0.0944-0.0956Columns46through54-0.0891-0.0905-0.0918-0.0931-0.0968-0.0980-0.0992Columns55through63-0.1004-0.1015-

9、0.1026-0.1037-0.1047-0.1057-0.1067-0.1077-0.1087Columns64through72-0.1096-0.1105-0.1114-0.1122-0.1130-0.1138-0.1146-0.1154-0.1161Columns73through81-0.1168-0.1175-0.1181-0.1187-0.1194-0.1199-0.1205-0.1210-0.1215Columns82through90-0.1220-0.1224-0.1229-0.1232-0.1236-0.1240-0.1243-0.1246-0.1249Columns91

10、through99-0.1251-0.1253-0.1255-0.1257-0.1259-0.1260-0.1261-0.1262-0.1262Columns100through101-0.1262-0.12621.0e-004*-0.0025-0.0037-0.0052-0.0218-0.0251-0.0286Columns1through9-0.0002-0.0007-0.00150.0069-0.0088-0.0109Columns10through18-0.0133-0.0159-0.01880.0323-0.0362-0.0403Columns19through27-0.0446-0

11、.0492-0.0539-0.0588-0.0639-0.0692-0.0747-0.0804-0.0863Columns28through36-0.0924-0.0986-0.1050-0.1116-0.1184-0.1253-0.1324-0.1397-0.1471Columns37through45-0.1547-0.1624-0.1703-0.1783-0.1865-0.1949-0.2034-0.2120-0.2208Columns46through54-0.2297-0.2388-0.2479-0.2572-0.2667-0.2762-0.2859-0.2957-0.3057Col

12、umns55through63- 0.3157-0.3258-0.3361-0.3782-0.3890-0.3998Columns64through72- 0.4108-0.4218-0.4330-0.4784-0.4899-0.5015Columns73through81- 0.5132-0.5249-0.5367-0.5846-0.5967-0.6088Columns82through90- 0.6210-0.6333-0.6456-0.6951-0.7075-0.7200Columns91through99- 0.7325-0.7451-0.7576-0.3465-0.3569-0.36

13、75-0.4442-0.4555-0.4669-0.5486-0.5606-0.5726-0.6579-0.6703-0.6827-0.7702-0.7828-0.7954-0.8080-0.8206-0.8332-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%绘弯矩图EditTextGoCellTo.olsDebugDezktopWindowHelpp31同£*：+ffleg-1.0+.11x

14、喊碰GtL-L=l;P=1000;Ll=l;蜓合出常数2-E=200*10*9；1=2*10*-5；3 -x=lizisp&ce(0L,101);dx=L/10口；鬻乂分1口口段4 -nl=Ll/dM+l;%确定k=LI,界对向的下标5 -N1=_F*(Ll-s(1;%第一段弯矩赋值6 -M2=zeros(L101-nlk%第二段弯矩呢:值(全为雪)7 -M=MLM2;抬全梁的弯矩S-A=curnwumCE*工)%对弯矩枳分求转角勺-Y-匚umwum(A)*dx%对转角枳分求挠度10 -subplotC3flDplot(je,grid找绘弯矩图11 -与ubpl”1,2)11PLota

15、,A)/grid%筵弯矩图12 -subplot0,L3),Plot(x,Y),grid抬绘弯矩图流程图n1=L"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)*dx结束(2)计算积分1一“解题思路：exp(-xA2)是不可积函数，matlab中int积分显示不出积分结果，用到vpa函数控制其结果精度，表示出来。程序：>>symsx>>t=vpa(int(exp(-x.A2)./(1+x.A2),-inf,+inf),5)结

16、果：t=1.3433J1工解题思路：先建立内置函数，然后调用quad函数求积分。程序：>>fun=(x)tan(x)./x.A(0.7);>>quad('fun',0,1)结果：ans=0.9063(普，L;解题思路：先建立内置函数，然后调用quad函数求积分。程序：>>fun=inline('exp(x)./(1-x.A2).A0.5');>>quad('fun',0,1)结果：ans=3.1044“(l+彳+/6dm为彳口+/<2x解题思路：这是一个二重积分，一般的方法是把二重积分化为二次积分，但由于二次积分命