《数学实验求微分方程的解》由会员分享，可在线阅读，更多相关《数学实验求微分方程的解（25页珍藏版）》请在金锄头文库上搜索。
1、实验四求微分方程的解1q 自牛顿发明微积分以来，微分方程在描述事物运动规律上已发挥了重要的作用。实际应用问题通过数学建模所得到的方 程，绝大多数是微分方程。q 由于实际应用的需要，人们必须求解微分方程。然而能够求得解析解的微分方程十分有限，绝大多数微分方程需要 利用数值方法来近似求解。q 本实验主要研究如何用 Matlab 来计算微分方程（组）的数值解，并重点介绍一个求解微分方程的基本数值解法 Euler折线法。问题背景和实验目的2q 考虑一维经典初值问题l 基本思想：用差商代替微商根据 Talyor 公式，y(x) 在点 xk 处有Euler 折线法3初值问题的Euler折线法q 具体步骤：
2、等距剖分：步长：l 分割求解区间分割求解区间，差商代替微商，解代数方程为分割点l 差商代替微商得迭代格式：k = 0, 1, 2, ., n-1yk 是 y (xk) 的近似4Euler 折线法举例例：用 Euler 法解初值问题取步长 h = (2 - 0)/n = 2/n，得差分方程当 h=0.4，即 n=5 时，Matlab 源程序见 fuluA.m解：5Euler 折线法源程序clear; f = inline(y+2*x/y2,x,y); a = 0; b = 2; n = 5; h = (b-a)/n; x = a : h : b; y(1) = 1; for i = 1 : ny
3、(i+1) = y(i) + h*f(x(i),y(i); end plot(x,y,rs-);6Euler折线法举例（续）解析解：解析解近似解y=1/3*(-18-54*x+45*exp(3*x)(1/3)7Runge-Kutta 方法q 为了减小误差，可采用以下方法：l 让步长 h 取得更小一些；l 改用具有较高精度的数值方法，如：q 龙格-库塔方法Runge-Kutta (龙格-库塔) 方法l 是一类求解常微分方程的数值方法l 有多种不同的迭代格式8Runge-Kutta 方法q 用得较多的是 四阶R-K方法（教材第 98 页）其中9四阶 R-K 方法源程序clear f = inlin
4、e(y+2*x/y2,x,y); a = 0; b = 2; n = 5; h = (b-a)/n; x = a : h : b; y(1) = 1; for i = 1 : nL1 = 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 点的近似值为 %fn,y(n+1); plot(x,y,rs-);10Runge-Kutta 方
5、法11Euler 法与 R-K法误差比较12Matlab 解初值问题函数q 用 Maltab自带函数 解初值问题l 求解析解：dsolvel 求数值解： ode45、ode23、ode113、ode23t、ode15s、ode23s、ode23tb13dsolve 求解析解q dsolve 的使用y=dsolve(eq1,eq2, . ,cond1,cond2, . ,v)其中 y 为输出的解， eq1、eq2、. 为微分方程，cond1、 cond2、. 为初值条件，v 为自变量。例 1：求微分方程 的通解，并验证。 y=dsolve(Dy+2*x*y=x*exp(-x2),x) syms
6、x; diff(y)+2*x*y - x*exp(-x2)14dsolve 的使用q 几点说明l 如果省略初值条件，则表示求通解；l 如果省略自变量，则默认自变量为 t dsolve(Dy=2*x,x); % dy/dx = 2x dsolve(Dy=2*x); % dy/dt = 2xl 若找不到解析解，则返回其积分形式。l 微分方程中用 D 表示对 自变量 的导数，如：Dy y； D2y y； D3y y15dsolve 举例例 2：求微分方程 在初值条件 下的特解，并画出解函数的图形。 y=dsolve(x*Dy+y- exp(x)=0,y(1)=2*exp(1),x) ezplot(y
7、);16dsolve 举例例3：求微分方程组 在初值条件 下的特解，并画出解函数的图形。x,y=dsolve(Dx+5*x+y=exp(t),Dy-x-3*y=0, . x(0)=1, y(0)=0, t) ezplot(x,y,0,1.3);注：解微分方程组时，如果所给的输出个数与方程个数相同 ，则方程组的解按词典顺序输出；如果只给一个输出，则输 出的是一个包含解的结构（structure）类型的数据。17dsolve 举例例：x,y=dsolve(Dx+5*x=0,Dy-3*y=0, .x(0)=1, y(0)=1,t)r = dsolve(Dx+5*x=0,Dy-3*y=0, .x(0)
8、=1, y(0)=1,t)这里返回的 r 是一个 结构类型 的数据r.x %查看解函数 x(t) r.y %查看解函数 y(t)只有很少一部分微分方程（组）能求出解析解。 大部分微分方程（组）只能利用数值方法求数值解。dsolve的输出个数只能为一个 或 与方程个数相等18数值求解T,Y = solver(odefun,tspan,y0)其中 y0 为初值条件，tspan为求解区间；Matlab在数值求解 时自动对求解区间进行分割，T (列向量) 中返回的是分割点 的值(自变量)，Y (数组) 中返回的是这些分割点上的近似解，其列数等于因变量的个数。 solver 为Matlab的ODE求解器
9、（可以是 ode45、ode23、 ode113、ode15s、ode23s、ode23t、ode23tb）没有一种算法可以有效地解决所有的 ODE 问题，因此 MATLAB 提供了多种ODE求解器，对于不同的ODE，可以 调用不同的求解器。19Matlab的ODE求解器求解器 ODE类 型特点说明ode45 非刚性 单步法；4，5 阶 R- K 方法；累计截断 误差为 (x)3大部分场合的首 选方法ode23 非刚性 单步法；2，3 阶 R- K 方法；累计截断 误差为 (x)3使用于精度较低 的情形ode11 3非刚性 多步法；Adams算法 ；高低精度均可到 10-310-6计算时间比
10、ode45 短ode23 t适度刚 性采用梯形算法适度刚性情形ode15 s刚性多步法；Gears 反 向数值微分；精度 中等若 ode45 失效 时，可尝试使 用 ode23 s刚性单步法；2 阶 Rosebrock 算法； 低精度当精度较低时， 计算时间比 ode15s 短 ode23 tb刚性梯形算法；低精度当精度较低时， 计算时间比 ode15s短20参数说明odefun 为显式常微分方程，可以用命令 inline 定义，或在函数文件中定义，然后通过函数句柄调用。fun=inline(-2*y+2*x2+2*x,x,y); x,y=ode23(fun,0,0.5,1);注：也可以在 t
11、span 中指定对求解区间的分割，如：x,y=ode23(fun,0:0.1:0.5,1); % x=0:0.1:0.5T,Y = solver(odefun,tspan,y0)求 的数值解，求解区间为 0,0.5例 ：21数值求解举例如果需求解的问题是高阶常微分方程，则需将其化为一阶常 微分方程组，此时必须用函数文件来定义该常微分方程组。令 求解 Ver der Pol 初值问题例 ：22数值求解举例l 先编写函数文件 verderpol.mfunction xprime=verderpol(t,x) global mu; xprime=x(2); mu*(1-x(1)2)*x(2) - x
12、(1);l 再编写脚本文件 vdpl.m，在命令窗口直接运行该文件 clear; global mu; mu=7; y0=1; 0; t,x=ode45(verderpol, 0,40, y0); % t,x=ode45(verderpol, 0,40, y0); plot(t,x(:,1);23求解微分方程小结q Matlab 函数l 求解析解（通解或特解），用 dsolve l 求数值解（特解），用 ode45、ode23 .q Matlab 编程l Euler 折线法 l Runga-Kutta 方法24l 教材 P97：练习 1 7上机作业要求写实验报告l 将所编写的程序分别命名为 hw41.m, hw42.m, hw43.m, hw44.m（主文件）,hw44fun.m（函数文件）, hw45.m, hw46.m（RK法解第5题）, hw47.m（第6、7题写在一起）l 将所有 M 文件作为附件，发给 mhjssystem.maill 邮件主题为：机号-学号-姓名；其中机号为 两位数l 三个字段之间用英文状态下的减号链接q 上机要求每个 M文件的第一行为：% 机号-学号-姓名25
