MATLAB实验四求微分方程的解.pptxVIP

MATLAB实验四求微分方程的解会计学问题背景和实验目的第1页/共24页 自牛顿发明微积分以来，微分方程在描述事物运动规律上已发挥了重要的作用。实际应用问题通过数学建模所得到的方程，绝大多数是微分方程。 由于实际应用的需要，人们必须求解微分方程。然而能够求得解析解的微分方程十分有限，绝大多数微分方程需要利用数值方法来近似求解。 本实验主要研究如何用 Matlab 来计算微分方程（组）的数值解，并重点介绍一个求解微分方程的基本数值解法－－Euler折线法。第2页/共24页Euler 折线法 考虑一维经典初值问题 基本思想：用差商代替微商根据 Talyor 公式，y(x) 在点 xk 处有 为分割点等距剖分：步长：得方程组：初值问题的Euler折线法第3页/共24页 具体步骤：分割求解区间，差商代替微商，解代数方程 分割求解区间 差商代替微商k = 0, 1, 2, ..., n-1yk 是 y (xk) 的近似第4页/共24页Euler 折线法举例例：用 Euler 法解初值问题解：取步长 h = (2 - 0)/n = 2/n，得差分方程当 ，即 n=5 时，Matlab 源程序见 Euler 折线法源程序第5页/共24页clearf=sym(y+2*x/y^2);a=0; b=2;h=0.4;n=(b-a)/h+1; % n=(b-a)/h;x=0; y=1;szj=[x,y];for i=1:n-1 % i=1:n y=y+h*subs(f,{x,y},{x,y}); x=x+h; szj=[szj;x,y];endszjplot(szj(:,1),szj(:,2),or-) Euler折线法举例（续）第6页/共24页解析解：y=1/3*(-18-54*x+45*exp(3*x))^(1/3)解析解近似解第7页/共24页Runge-Kutta 方法 为了减小误差，可采用以下方法： 让步长 h 取得更小一些； 改用具有较高精度的数值方法：Runge-Kutta (龙格-库塔) 方法 龙格-库塔方法 是一类求解常微分方程的数值方法 有多种不同的迭代格式其中第8页/共24页Runge-Kutta 方法 用得较多的是 四阶R-K方法（教材第 98 页）四阶 R-K 方法源程序第9页/共24页clear;f=sym(y+2*x/y^2);a=0; b=2; h=0.4;n=(b-a)/h+1; % n=(b-a)/h;x=0; y=1; szj=[x,y];for i=1:n-1 % i=1:n l1=subs(f,{x,y},{x,y}); l2=subs(f,{x,y},{x+h/2,y+l1*h/2}); l3=subs(f,{x,y},{x+h/2,y+l2*h/2}); l4=subs(f,{x,y},{x+h,y+l3*h}); y=y+h*(l1+2*l2+2*l3+l4)/6; x=x+h; szj=[szj;x,y];endplot(szj(:,1),szj(:,2), dg-)第10页/共24页Runge-Kutta 方法Euler 法与 R-K法误差比较第11页/共24页Matlab 解初值问题第12页/共24页 用 Maltab自带函数 解初值问题 求解析解：dsolve 求数值解： ode45、ode23、 ode113、ode23t、ode15s、 ode23s、ode23tb第13页/共24页dsolve 求解析解 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)第14页/共24页dsolve 的使用 几点说明 微分方程中用 D 表示对 自变量 的导数，如：Dy y； D2y y； D3y y 如果省略初值条件，则表示求通解； 如果省略自变量，则默认自变量为 t dsolve(Dy=2*x,x); ％ dy/dx = 2xdsolve(Dy=2*x);％ dy/dt = 2x 若找不到解析解，则返回其积分形式。第15页/共24页dsolve 举例例 2：求微分方程 在初值条件 下的特解，并画出解函数的图形。 y=dsolve(x*Dy+y-exp(x)=0,y(1)=2*exp(1),x)