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

欧拉近似方法求常微分方程 朱翼 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;end if nargin6,trace=0;end t=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,y end while (ttfinal)(t+ht) %主循环 if t+htfinal,h=tfinal-t;end % Compute the slopes f=feval(ypfun,t,y);f=f(:); %估计误差并设定可接受误差 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 klength(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 size if delta~=0.0 h=min(hmax,0.9*h*(tau/delta)^pow); end end if (ttfinal) dish(Singularity likely.) t end tout=tout(1:k); yout=yout(1:k,:); 流程图： 举例： 用欧拉法求y’=-y+x+1，y（0）=1。 解题思路： 首先建立例子所给函数的函数文件，调用函数myeuler，利用程序求解方程。将欧拉法解和精确解比较，由方程f=-y+x+1可得到其精确解为y（x）=x+exp(-x)。最后在同一图窗中分别画出两图。 程序代码： f．m function 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(欧拉法,精确解) 例题流程图： 2、编程解决以下科学计算问题 （1） 解题思路：◆建模： 材料力学中从弯矩求转角要经过一次不定积分, 而从转角求挠度又要经过一次不定积分, 在MATLAB中这却是非常简单的问题.可用cumsum函数作近似的不定积分,还可用更精确的函数cumtrapz来做不定积分。本题用cumsum函数来做.解题的关键还是在于正确地列写弯矩方程。本例中弯矩为 程序代码： L=1; P=1000; L1=1; %给出常数 E = 200*10^9; 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 % 对转角积分求挠