资源预览内容
第1页 / 共12页
第2页 / 共12页
第3页 / 共12页
第4页 / 共12页
第5页 / 共12页
第6页 / 共12页
第7页 / 共12页
第8页 / 共12页
第9页 / 共12页
第10页 / 共12页
亲,该文档总共12页,到这儿已超出免费预览范围,如果喜欢就下载吧!
资源描述
第48讲 微分方程数值解,数学建模,描述自然世界客观规律的模型很多都是微分方程模型,但是很多都没有解析解,需要求出其数值解以满足需求。,假设其存在唯一且足够光滑的解,所谓数值解,就是寻求解在一些离散点,也被称之为节点, 处的近似值 。相邻2个节点的间距 ,一般取为常数 , 此时 。,(1),数学建模,Euler法,首先对(1)中的微分方程两端在区间 上积分,得到 利用左矩形公式计算右端积分,再用近似值 代替 ,则,Euler方法,数学建模,Euler法几何意义,为方程的精确解 ,在区间 上,用过点 、以 为斜率的直线 ,近似代替 ,用该直线与直线 的交点的纵坐标近似代替 。,Euler法也称为折线法,数学建模,例:求解,表11.3 用Euler方法计算得到的结果,称方法具有一阶精度。,数学建模,利用右矩形公式计算(2)右端的积分,再用近似值 代替 ,则,向后Euler公式,向后Euler公式是一个隐式公式,无法直接计算出 ,所以大多数时候,应用该公式要结合Euler公式计算,先用Euler公式计算出中间值(预估值)。,再用向后Euler公式计算校正值,如果 收敛,则取 。向后Euler公式仍然只具有一阶精度。,数学建模,利用Matlab求解微分方程, syms y mu; y=dsolve(D2y+mu*(y2-1)*Dy+y=0) Warning; Compact, analytic solution could not be found. dsolve命令不能够求出方程的解析解,dsolve 求解常微分方程的解析解 语法: dsolve(f1,f2,.,fm),f1,f2,.,fm为微分方程,也可以是初值条件。Matlab中采用Dky表示y对x的k阶导数。,例: 求解 y=dsolve(D2y=sin(t) y= sin(t)+C1*t+C2 这里C1和C2为任意常数,例: 求解,数学建模,利用Matlab求解微分方程,ode45/ode23 求解常微分方程数值的函数 语法: t,x=ode45ode23(Fun,t_0,t_f,x_0),Fun表示指定的M文件给出的常微分方程函数;t_0,t_f为自变量取 值;x_0表示方程的初始条件。Ode45和ode23分别用4(5)阶和2(3)阶Runge-Kutta方法求解。,例: 求解,建立函数文件fun2.m function xdot=fun2(t,theta) g=9.8; l=1; m=1; lamda=0.1; xdot=theta(2);-g*sin(theta(1)/1-lamda*theta(2)/m;,数学建模,Matlab求解微分方程,command window中键入: theta0=pi/12,0; t,theta=ode23(fun2,0,100,theta0); plot(t,theta(:,1),这个方程组是单摆在粘性介质中摆动的模型,可以看出单摆的运动最后趋于静止。,数学建模,Matlab求解微分方程,例: 求解 其中 。,首先编写M文件lorenzeq.m来描述方程组。 function xdot=lorenzeq(t,x) xdot=-8/3*x(1)+x(2)*x(3);-10*x(2)+10*x(3);-x(1)*x(2)+28*x(2)-x(3); 然后在命令窗口设定求解区域并给定初值进行求解 t_final=100;x0=0;0;1e-10; t,x=ode45(lorenzeq,0,t_final,x0); plot(t,x);,数学建模,Thank you,微分方程数值解,
网站客服QQ:2055934822
金锄头文库版权所有
经营许可证:蜀ICP备13022795号 | 川公网安备 51140202000112号