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

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 nargin5,tol=1.e-3; endif nar

2、gin6,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 traceclc,t,h,yendwhile (tt)if t+htfinal,h=tfinal-t; end % Compute the slopes f=feval(ypfun,t,y);f=f(:);%估计误差并设定可接受误差 delta=norm(h*f, inf ); tau=tol

3、*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 trace home,t,h,y end% Update the step sizeif delta=0.0 h=min(hmax,0.9*h*(tau/delta)pow); endendif (ttfinal)dish( Singularity likely. )tend tout=tout(1:k)

4、; yout=yout(1:k,:);流程图：开始Po3YNNtrace=NNttfinalNYNklengt(tout)f)1.0)t+ht Ynodrmel(tya, in=ft)a1u.tout,yout=myeuler(ypfun,t0,tfinal,y0,tol,trace)Y t=t0;hmax=(tfinal-t)/16;h=hmax/8;y=y0(:);chunk=128;nargin6tout=zeros(chunk,1);yout=zeros(chunk,length(y);k=1;tout(k)=t;yout(k,:)=y.;narginx,y=myeuler(f,0,1

5、,1); %利用程序求解方程 y1=x+exp(-x);%方程 f=-y+x+1 的精确解plot(x,y,-b,x,y1,-r) %在同一图窗将欧拉法解和精确解的图画出 legend(欧 拉法 ,精确解 )例题流程图： 开始f=f(x,y)f=-y+x+1;调用函数myeuler ，x,y=myeuler( f ,0,1,1);求出y结1果=x+exp(-x)2、1)解题思路 ：建模： 材料力学中从弯矩求转角要经过一次不定积分 , 而从转角求挠度又要经过一次不定积分 , 在 MATLAB中这却 是非常简单的问题 .可用 cumsum函数作近似的不定积分 ,还可 用更精确的函数 cumtrap

6、z 来做不定积分。本题用 cumsum 函 数来做.解题的关键还是在于正确地列写弯矩方程。本例中弯 矩为程序代码： 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 = cums

7、um(A)*dx% 对转角积分求挠度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.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 t

8、hrough 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.08350.0849 -0.0863 -0.0877Columns 46 through 54-0.0891 -0.0905 -0.0918-0.0931-0.0944-0.09560.0968 -0.0980 -0.0992Columns 55 through 63-0.1004 -0.1015 -0.1026-0.103

9、7-0.1047-0.10570.1067 -0.1077 -0.1087Columns 64 through 72-0.1096 -0.1105 -0.1114-0.1122-0.1130-0.11380.1146 -0.1154 -0.1161Columns 73 through 81-0.1168 -0.1175 -0.1181-0.1187-0.1194-0.11990.1205 -0.1210 -0.1215Columns 82 through 90-0.1220 -0.1224 -0.1229-0.1232-0.1236-0.12400.1243 -0.1246 -0.1249Co

10、lumns 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.12621.0e-004 *Columns 1 through 9-0.0002 -0.0007 -0.0015-0.0025-0.0037-0.00520.0069 -0.0088 -0.0109Columns 10 through 18-0.0133 -0.0159 -0.0188-0.0218-0.0251-0.02860.0323 -0.0362 -

11、0.0403Columns 19 through 27-0.0446 -0.0492 -0.0539-0.0588-0.0639-0.06920.0747 -0.0804 -0.0863Columns 28 through 36-0.0924 -0.0986 -0.1050-0.1116-0.1184-0.12530.1324 -0.1397 -0.1471Columns 37 through 45-0.1547 -0.1624 -0.1703-0.1783-0.1865-0.19490.2034 -0.2120 -0.2208Columns 46 through 54-0.2297 -0.2

12、388 -0.2479-0.2572-0.2667-0.27620.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.5015Columns 73 through 81-0.5132 -0.5249 -0.5367-0.5486-0.5606-0.5726-0.58

13、46 -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 101-0.8459 -0.8585 subplot(3,1,1),plot(x,M),grid % 绘弯矩图 subplot(3,1,2),plot(x,A),grid%

14、 绘弯矩图subplot(3,1,3),plot(x,Y),grid% 绘弯矩图流程图开始解题思路 ：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.2,y,-sq