数学实验：求微分方程及解.ppt

实验四 问题背景和实验目的 Euler 折线法 初值问题的Euler折线法 Euler 折线法举例 Euler 折线法源程序 Euler折线法举例（续） Runge-Kutta 方法 Runge-Kutta 方法 四阶 R-K 方法源程序 Runge-Kutta 方法 Euler 法与 R-K法误差比较 Matlab 解初值问题函数 dsolve 求解析解 dsolve 的使用 dsolve 举例 dsolve 举例 dsolve 举例 数值求解 Matlab的ODE求解器 参数说明 数值求解举例 数值求解举例 求解微分方程小结 上机作业 * * 求微分方程的解 自牛顿发明微积分以来，微分方程在描述事物运动规律上已发挥了重要的作用。实际应用问题通过数学建模所得到的方程，绝大多数是微分方程。 由于实际应用的需要，人们必须求解微分方程。然而能够求得解析解的微分方程十分有限，绝大多数微分方程需要利用数值方法来近似求解。 本实验主要研究如何用 Matlab 来计算微分方程（组）的数值解，并重点介绍一个求解微分方程的基本数值解法－－Euler折线法。 考虑一维经典初值问题 基本思想：用差商代替微商 根据 Talyor 公式，y(x) 在点 xk 处有 具体步骤： 等距剖分： 步长： 分割求解区间 分割求解区间，差商代替微商，解代数方程 为分割点 差商代替微商 得迭代格式： k = 0, 1, 2, ..., n-1 yk 是 y (xk) 的近似 例：用 Euler 法解初值问题 取步长 h = (2 - 0)/n = 2/n，得差分方程 当 h=0.4，即 n=5 时，Matlab 源程序见 fuluA.m 解： clear; f = inline(y+2*x/y^2,x,y); a = 0; b = 2; n = 5; h = (b-a)/n; x = a : h : b; y(1) = 1; for i = 1 : n y(i+1) = y(i) + h*f(x(i),y(i)); end plot(x,y,rs-); 解析解： 解析解 近似解 y=1/3*(-18-54*x+45*exp(3*x))^(1/3) 为了减小误差，可采用以下方法： 让步长 h 取得更小一些； 改用具有较高精度的数值方法，如： 龙格-库塔方法 Runge-Kutta (龙格-库塔) 方法 是一类求解常微分方程的数值方法 有多种不同的迭代格式 用得较多的是 四阶R-K方法（教材第 98 页） 其中 clear f = inline(y+2*x/y^2,x,y); a = 0; b = 2; n = 5; h = (b-a)/n; x = a : h : b; y(1) = 1; for i = 1 : n L1 = f(x(i), y(i)); L2 = f(x(i)+h/2, y(i)+L1*h/2); L3 = f(x(i)+h/2, y(i)+L2*h/2); L4 = f(x(i)+h, y(i)+L3*h); y(i+1) = y(i) + h*(L1+2*L2+2*L3+L4)/6; end fprintf(y 在 x=2 点的近似值为 %f\n,y(n+1)); plot(x,y,rs-); 用 Maltab自带函数 解初值问题 求解析解：dsolve 求数值解：  ode45、ode23、 ode113、ode23t、ode15s、 ode23s、ode23tb dsolve 的使用 y=dsolve(eq1,eq2, ... ,cond1,cond2, ... ,v) 其中 y 为输出的解， eq1、eq2、... 为微分方程，cond1、cond2、... 为初值条件，v 为自变量。 例 1：求微分方程 的通解，并验证。 y=dsolve(Dy+2*x*y=x*exp(-x^2),x) syms x; diff(y)+2*x*y - x*exp(-x^2) 几点说明 如果省略初值条件，则表示求通解； 如果省略自变量，则默认自变量为 t dsolve(Dy=2*x,x); % dy/dx = 2x dsolve(Dy=2*x); % dy/dt = 2x 若找不到解析解，则返回其积分形式。 微分方程中用 D 表示对 自变量 的导数，如： Dy y； D2y y； D3y y 例 2：求微分方程 在初值条件