19章 微分方程数值解.docVIP

第19章 常微分方程初值问题数值解 Euler方法 Euler方法算法 输入函数f(x,y)、初值y0、自变量区间端点a,b步长h 计算节点数n和节点x k 用Euler公式y n+1= yn+h f(xn,yn) 求数值解 Euler方法程序 Clear[x,y,f] f[x_,y_]= Input[函数f(x,y)=] y0=Input[初值y0 =] a=Input[左端点a=] b=Input[右端点b=] h=Input[步长h=] n=(b-a)/h; For[i=0,in,i++, xk=a+i*h; y1=y0+h*f[xk,y0]; Print[y(,xk+h//N,)=,y1//N]; y0=y1] 说明：本程序用Euler公式求常微方程初值问题数值解。程序执行后，按要求通过键盘依次输入输入函数f(x,y)、初值y0、自变量区间端点a,b步长h后，计算机则给出常微方程初值问题数值解。 程序中变量说明 f[x,y]: 存放函数f(x,y) y0: 存放初值y0及数值解 a：存放自变量区间左端点 b: 存放自变量区间右端点 n: 存放节点个数 h: 存放节点步长 xk：存放节点xi y1: 存放数值解 注：（1）语句Print[y(,xk+h//N,)=,y1//N]是将数值解用6位有效数字显示出来，如果要显示n位数有效数字，可以将语句改为Print[y(,xk+h//N,)=,N[y1,n]] （2）Mathematica中有求微分方程初值问题数值解的命令，形式为: NDSolve[ {y’[x]==f[x,y[x]],y[a]==y0}, y, {x, a, b}] 由NDSolve命令得到的解是以{{未知函数名-InterpolatingFunction[range, ]}}的形式给出的,其中的InterpolatingFunction[range, ]是所求的插值函数表示的数值解, range就是所求数值解的自变量范围。 例题与实验 用Euler方法求初值问题 的数值解。取步长h=0.1,并在一个坐标系中画出数值解与准确解的图形。 解：执行Euler方法程序后，在输入的窗口中按提示分别输入y-2x/y、1、0、1、0.1，每次输入后用鼠标点击窗口的“OK”按扭，计算机在屏幕上给出如下数值解结果。 y(0.1)=1.1 y(0.2)=1.19182 y(0.3)=1.27744 y(0.4)=1.35821 y(0.5)=1.43513 y(0.6)=1.50897 y(0.7)=1.58034 y(0.8)=1.64978 y(0.9)=1.71778 y(1.)=1.78477 本题的准确解为y(x)= (1 + 2 x)1/2，在一个坐标系中画出数值解与准确解的图形如下。 图中点图是数值解，曲线为准确解。从图中可以知道数值解与准确解比较接近，但还是有误差的。 改进的Euler方法 改进的Euler方法算法 1. 输入函数f(x,y)、初值y0、自变量区间端点a,b步长h 2.计算节点数n和节点x k 3.用改进的Euler公式求数值解 改进的Euler方法程序 Clear[x,y,f] f[x_,y_]= Input[函数f(x,y)=] y0=Input[初值y0 =] a=Input[左端点a=] b=Input[右端点b=] h=Input[步长h=] n=(b-a)/h; For[i=0,in,i++, xk=a+i*h; y1=y0+h*f[xk,y0]; y1=y0+h/2*(f[xk,y0]+f[xk+h,y1]); Print[y(,xk+h//N,)=,y1//N]; y0=y1] 注：语句h=Input[步长h=]的步长输入应该用带小数点的数输入，这样可以加快计算速度，否则，由于Mathematica软件本身的原因，可能会出现计算时间过长等计算不出结果的情况。 说明：本程序用改进的Euler公式求常微方程初值问题数值解。程序执行后，按要求通过键盘依次输入输入函数f(x,y)、初值y0、自变量区间端点a,b步长h后，计算机则给出常微方程初值问题数值解。 程序中变量说明： f[x,y]: 存放函数f(x,y) y0: 存放初值y0及数值解 a：存放自变量区间左端点 b: 存放自变量区间右端点 n: 存放节点个数 h: 存放节点步长 xk：存放节点xi y1: 存放数值解 例题与实验 用改进的Euler方法求初值问题 的数值解。取步长h=0.1和0.01计算，并与准确解比较。 解：执行改进的Euler方法程序后，在输入的窗口中按提示分别输入-2x*y^2、1、0、1、0.1，每次输入后用鼠标点击窗口的“OK”按扭，计算机在屏幕上给出如下数值解结果