【作者声明】


  本文所有文字均为作者原创,所有图片均为作者本人真实拍摄或制作。

  版权所有,仅供阅读欣赏,严禁任何单位或个人以任何形式转载、复制、引用、抄袭、截图、模仿、翻译本评测的部分或全部内容(包括但不限于文字、图片)。

  作者保留所有权利。

  请尊重作者劳动成果,谢谢合作。


前言

  一阶微分方程是最常见的一类微分方程。我们在《高等数学》等课程中学过一些微分方程的解法,例如齐次方程、线性方程等等。然而在实际应用中,我们更关心的是微分方程在某个初始条件下的特解,并且很多微分方程无法使用通常的初等积分法求解,这时运用数值法或几何法求解微分方程的优势就体现出来了。
这里运用几何法来求解微分方程,即先绘制出微分方程的斜率场,然后在斜率场中从初始条件开始绘制微分方程特解对应的积分曲线。其实,有很少的一部分图形计算器(一般是具有CAS功能的图形计算器)拥有内置的绘制微分方程斜率场的功能,然而大多数图形计算器却不具备这种功能,因此本文利用图形计算器的编程绘图功能,编写程序绘制微分方程斜率场,进而再画出积分曲线。

数学原理

一、微分方程斜率场

  当一阶微分方程y’=f(x,y)不能或无法轻易地运用一般方法求出解的时候,可以采用近似解法求解。设该方程右边的函数f(x,y)在某区域D内有定义,那么过D内的每一点M(x,y)都可以作一条以f(x,y)为斜率的直线,即对于D内的每一点,都确定一个斜率值与之对应,即该微分方程在区域D内确定了一个斜率场。以点M为起点的或点M附近的一个线段,能直观地表示D内过该点的积分曲线的斜率。如果D内的一条曲线在任意一点的切线斜率都和斜率场在该点处的线段斜率一致,那么这条曲线就是该微分方程的一条积分曲线。
  在这篇文章中,我们不涉及关于一阶微分方程初值问题是否有解、解是否唯一的讨论,总是假定所讨论的一阶微分方程初值问题的解都是存在且唯一的。

二、常微分方程初值问题的数值解法

  常微分方程初值问题的数值解法有很多,例如欧拉(Euler)法、龙格-库塔(Runge-Kutta)方法、亚当斯(Adams)方法等等,本文使用四阶龙格-库塔方法作为积分曲线绘制的算法。算法如下:

对于微分方程y’=f(x,y),有初始条件y(x0)=y0,那么:y[n+1]=y[n]+h(k1+2·k2+2·k3+k4)/6,其中
k1=f(x[n],y[n])
k2=f(x[n]+h/2,y[n]+h·k1/2)
k3=f(x[n]+h/2,y[n]+h·k2/2)
k4=f(x[n]+h,y[n]+h·k3)
h为迭代的步长值。

程序方案

一、程序文件

  为了保护程序代码不被修改,我们将程序分为绘图主程序(DIFFGRPH)与微分方程程序(DIFFEQN)两部分,绘图主程序用于绘制斜率场及积分曲线,作为程序入口调用,无需修改;微分方程程序仅需要写入绘制的微分方程,根据需要进行修改。

二、斜率场绘制

  按照一般图形计算器的屏幕大小,将屏幕视窗区域分为8行,然后按照行宽划分纵列,这样屏幕就被划分为很多个小正方形区域,以小正方形区域中心为起点,绘制过该点的斜率场线段。

三、绘制给定初始条件的积分曲线

  在斜率场绘制完成后,输入初始条件,迭代的步长值默认使用计算器的绘制间隔,即可开始积分曲线的绘制。为了突出显示初始条件,可以以初始条件点为圆心,画一个半径很小的圆,然后向X值增大的方向绘制积分曲线,再回到初始条件,向X值减小的方向绘制积分曲线,使得曲线贯穿整个屏幕。

四、颜色分配

  如果是彩屏图形计算器,最好在绘图的时候加上颜色指令,让斜率场线段、初始条件、积分曲线分别以三种颜色显示,得到更好的绘图效果。

程序算法

  我们分别以CASIO的fx-CG50图形计算器和德州仪器的TI-84 Plus为例进行说明,求解的微分方程为y’=sin x cos y,初始条件为x=-π/2,y=0。

CASIO fx-CG50

绘图主程序(DIFFGRPH):
 Cls
 Rad
 (Xmax-Xmin)÷32→E
 (Ymax-Ymin)÷16→I
 For Xmin+E→A To Xmax-E Step 2E
 For Ymax-I→B To Ymin+I Step -2I(使用For循环划分正方形小区域)
 A→X
 B→Y
 Prog "DIFFEQN"
 Green F-Line A-0.5Icos tan⁻¹ P,B-0.5Isin tan⁻¹ P,A+0.5Icos tan⁻¹ P,B+0.5Isin tan⁻¹ P
 (使用绿色线以正方形小区域中心为基准,以该点为中心,按照该点的斜率值绘制一段线段)
 Next
 Next
 "X0="?→W
 "Y'(X0)="?→Z(输入初始条件)
 Z→R
 (Xmax-Xmin)/128→H
 Blue Circle W,Z,2H(用蓝色突出显示初始条件点)
 For W→A To Xmax Step H(向X值增大的方向绘制积分曲线)
 A→U
 Z→V
 A→X
 Z→Y
 Prog "DIFFEQN"
 P→J
 A+H÷2→X
 Z+HJ÷2→Y
 Prog "DIFFEQN"
 P→K
 Z+HK÷2→Y
 Prog "DIFFEQN"
 P→L
 A+H→X
 Z+HL→Y
 Prog "DIFFEQN"
 P→M
 Z+H(J+2K+2L+M)÷6→Z
 Red F-Line U,V,X,Z(用红色绘制积分曲线)
 Next
 R→Z
 For W→A To Xmin Step -H(向X值减小的方向绘制另一半积分曲线)
 A→U
 Z→V
 A→X
 Z→Y
 Prog "DIFFEQN"
 P→J
 A-H÷2→X
 Z-HJ÷2→Y
 Prog "DIFFEQN"
 P→K
 Z-HK÷2→Y
 Prog "DIFFEQN"
 P→L
 A-H→X
 Z-HL→Y
 Prog "DIFFEQN"
 P→M
 Z-H(J+2K+2L+M)÷6→Z
 Red F-Line U,V,X,Z(用红色绘制积分曲线)
 Next

微分方程程序(DIFFEQN)
 sin Xcos Y→P
 (将f(x,y)表达式的值赋值给变量P)

TI-84 Plus

绘图主程序(DIFFGRPH):
  ClrDraw
  Radian
  (Xmax-Xmin)/24→E
  (Ymax-Ymin)/16→I
  For(A,Xmin+E,Xmax-E,2E)
  For(B,Ymax-I,Ymin+I,­2I)(使用For循环划分正方形小区域)
  A→X
  B→Y
  prgmDIFFEQN
  Line(A-0.5Icos(tan⁻¹(P)),B-0.5Isin(tan⁻¹(P)),A+0.5Icos(tan⁻¹(P)),B+0.5Isin(tan⁻¹(P)))
  (以正方形小区域中心为基准,以该点为中心,按照该点的斜率值绘制一段线段)
  End
  End
  Disp "X0="
  Input W
  Disp "Y'(X0)="
  Input Z(输入初始条件)
  (Xmax-Xmin)/96→H
  Circle(W,Z,2H)(突出显示初始条件点)
  Z→R
  For(A,W,Xmax,H)(向X值增大的方向绘制积分曲线)
  A→U
  Z→V
  A→X
  Z→Y
  prgmDIFFEQN
  P→J
  A+H/2→X
  Z+HJ/2→Y
  prgmDIFFEQN
  P→K
  Z+HK/2→Y
  prgmDIFFEQN
  P→L
  A+H→X
  Z+HL→Y
  prgmDIFFEQN
  P→M
  Z+H(J+2K+2L+M)/6→Z
  Line(U,V,X,Z)(绘制积分曲线)
  End
  R→Z
  For(A,W,Xmin,­H)(向X值减小的方向绘制另一半积分曲线)
  A→U
  Z→V
  A→X
  Z→Y
  prgmDIFFEQN
  P→J
  A-H/2→X
  Z-HJ/2→Y
  prgmDIFFEQN
  P→K
  Z-HK/2→Y
  prgmDIFFEQN
  P→L
  A-H→X
  Z-HL→Y
  prgmDIFFEQN
  P→M
  Z-H(J+2K+2L+M)/6→Z
  Line(X,Z,U,V)(绘制积分曲线)
  End

微分方程程序(DIFFEQN):
  sin(X)cos(Y)→P
  (将f(x,y)表达式的值赋值给变量P)

程序运行

  本文提供的程序没有涉及到视窗的设置,运行程序前,应当先检查视窗范围,手动设置视窗到合适的范围。通常情况下,使用默认的初始窗即可。整幅图像绘制完成后,计算器会停留在图形界面,如果需要退出图形界面,按下[AC/ON]或[CLEAR]键即可。

程序运行效果:

CASIO fx-CG50:视窗为默认的初始窗

TI-84 Plus:视窗为ZDecimal

整段程序运行下来需要的时间并不多,基本上都在1分钟左右,比较适合实际应用。

总结

  我们又一次巧妙地利用图形计算器的编程功能,解决了微分方程斜率场绘图的问题。本文的程序中使用四阶龙格-库塔方法计算并绘制积分曲线的这一段程序,也可以稍作删减修改之后用于fx-5800P等无绘图功能的可编程计算器上,作为常微分方程初值问题的数值求解程序使用。从这一实例可以看到,图形计算器编程功能中自带的绘图命令(点、线、圆、法线、切线等等),与其他程序命令结合之后,能够很方便地完成复杂的图形绘制,从而解决相对复杂的数学问题。