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

欧拉近似方法求常微分方程 朱翼 1、编程实现以下科学计算算法，并举一例应用之。 “欧拉近似方法求常微分方程” 算法说明: 欧拉法是简单有效的常微分方程数值解法，欧拉法有多种形式的 算法，其中简单欧拉法是一种单步递推算法。其基本原理为 对简单的一阶方程的初值问题: ygy) 其中 y(x0)=y0 欧拉法等同于将函数微分转换为数值微分，rti欧拉公式可得 yn+i =yn+hf(xn/yn) 程序代码： function [tout, yout]=myeuler(ypfun, tO, tfinal, yO, tol, trace) %初始化 pow二1/3; if n arg in〈5, tol 二1. e~3; e nd if nargin〈6, trace二0;end t=tO; hmax=(tfinal~t)/16; h=hma.x/8; y=yO (:); chunk二128; tout=zeros(chunk, 1); yout=zeros(chunk, length(y)); k=l； tout(k)二t; yout (k,:)二y.; %绘图%主循环if trace %绘图 %主循环 clc, t, h, y end whi 1 e nal)(t+ht) if t+htfinal, h=tfinal~t;end % Compute the slopes f=fcval(ypfun, t, y);f=f(:); %估计误差并设定可接受误差 delta二norm(h*f, inf); tau=tol*max(norm(y,? inf，), 1. %当误差可接受时重写解 if delta=tau t二t+h; y二y+h*f； k=k+l; 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 si ze if delta^=0. 0 h=min(hmax, 0. 9*h*(tau/della) pow); encl end if (ttfinal) dish(J Singularity likely.) t end tout二tout(1:k); yout=yout(1:k,:); 流程图: t+htfi nal Y Y 举例： 用欧拉法求yJ-y+x+1, y (0) =1。 解题思路： 首先建立例子所给函数的函数文件，调用函数myeuler,利用程序 求解方程。将欧拉法解和精确解比较，由方程f二-y+x+l可得到其精确 解为y (x) =x+exp(-x)o最后在同一图窗屮分别画出两图。 程序代码: f. m function f二f(x, y) f二-y+x+l; ?[x,y]=myeuler(,f,,0,1,1); % 利用程序求解方程 ? yl=x+exp(-x); %方程 f二-y+x+l 的精确解 ?plot(x,y；-b,,x,yl；-r,) %在同一图窗将欧拉法解和精确解的图画出 ? legendf*欧拉法T精确解〕 例题流程图:  结束 2.编程解决以下科学计算问题 (1) Lfl 7-1-2：长为 L 的左端因 孰在离固定端5处施加力人求它的转角和挠度.设梁 E^ZOOXIO*N/m^l-SX^m4 为巳虬 解题思路：?建模：材料力学中从弯矩求转角要经过一次不定积分, 而从转角求挠度又要经过一次不定积分，在MATLAB中这却 是非常简单的问题?可用cumsum函数作近似的不定积分，还可 用更精确的函数cumtrapz来做不定积分。本题用cumsum函 数来做?解题的关键还是在于止确地列写弯矩方程。本例中弯 矩为(0xZJ (厶 xL) 矩为 (0xZJ (厶 xL) 转角“j;士此挠度心匚血 程序代码： ? L=l; P=1000; Ll=l;%给出常数 E = 200*10A9; l=2*10A-5; x = linspace(0,L,101); dx=L/100;%将 x 分 100 段 nl=Ll/dx+l; %确定x二LI处对应的下标 Ml = -P*( Ll-x(l:nl)); % 第一段弯矩赋值 M2 = zeros(l,101-nl); %第二段弯矩赋值(全为零) M = [M1,M2]; %全梁的弯矩 A = cumsum(M)*dx/(E*l) %对弯矩积分求转角 Y = cumsum(A)*dx %对转角积分求挠度 A = 1.0e-003 * Columns 1 through 9 -0.0146