MATLAB解微分方程

nullnullMATLAB ODE初值问题的数值解 PDE问题的数值解null问题提出 倒葫芦形状容器壁上的刻度问题.对于如图所示圆柱形状容器壁上的容积刻度,可以利用圆柱体体积公式其中直径D为常数. 对于几何形状不是规则的容器,比如倒葫芦形状容器壁上如何标出刻度呢?下 表 关于同志近三年现实表现材料材料类招标技术评分表图表与交易pdf视力表打印pdf用图表说话 pdf 是经过测量得到部分容器高度与直径的关系.nullx1o根据上表的数据,可以拟合出倒葫芦形状容器的图,建立如图所示的坐标轴,问题为如何根据任意高度x标出容器体积V的刻度,由微元思想 分析 定性数据统计分析pdf销售业绩分析模板建筑结构震害分析销售进度分析表京东商城竞争战略分析 可知H 0 0.2 0.4 0.6 0.8 1.0 D 0.04 0.11 0.26 0.56 1.04 1.17null只要求解上述方程,就可求出体积V与高度x之间的函数关系,从而可标出容器壁上容积的刻度,但问题是函数D(x)无解析表达式,无法求出其解析解. 因此, 得到如下微分方程初值问题null包含自变量、未知函数及未知函数的导数或微分的方程称为微分方程。 在微分方程中, 自变量的个数只有一个, 称为常微分方程。 自变量的个数为两个或两个以上的微分方程叫偏微分方程。 微分方程中出现的未知函数最高阶导数的阶数称为微分方程的阶数。 微分方程分类null常微分方程：(1)，(2)式称为初值问题.在实际应用中还经常需要求解常微分方程组：(3)式称为边值问题。null但能求解析解的常微分方程是有限的，大多数的常微分方程是给不出解析解的. 这个一阶微分方程就不能用初等函数及其积分来表达它的解。 例 例 nullnull事实上，从实际问题当中抽象出来的微分方程，通常主要依靠数值解法来解决。null微分方程数值 方法 快递客服问题件处理详细方法山木方法pdf计算方法pdf华与华方法下载八字理论方法下载 的基本思想 对常微分方程初值问题的数值解法，就是要算出精确解y(x)在区间a,b上的一系列离散节点 处的函数值 相邻两个节点的间距 称为步长，步 长可以相等，也可以不等。假定h为定数，称为定步长，这时节点可表示为null数值解法需要把连续性的问题加以离散化，从而求出离散节点的数值解。离散化 对常微分方程数值解法的基本出发点就是离散化。 其数值解法的基本特点：采用“步进式”， 即求解过程顺着节点排列的次序一步一步地向前推进，null描述这类算法，要求给出用已知信息 计算 的递推公式。 建立这类递推公式的基本方法是在这些节点上用数值积分、数值微分、泰勒展开等离散化方法，对初值问题 数值解和精确解数值解和精确解用数值方法求解初值问题，不是求出它的解析解或其近似解析式，而是给出它的解在某些离散节点上的近似值。用y(x)表示问题的准确解y(x0), y(x1),y(xN) 表示解y(x)在节点x0, x1,…, xN处的准确值 y0,y1,…,y N表示数值解，即问题的解y(x) 在相应节点处的近似值。单步法和多步法单步法和多步法单步法：在计算yi+1 时只利用y i 多步法：在计算yi+1 时不仅利用y i , 还要利用 yi−1, yi−2,…, k步法：在计算yi+1 时要用到yi,yi−1,…,yi−k+1显式格式可写成:yk+1=yk+hΦf(xk，yk;h） 隐式格式:yk+1=yk+hΦf（xk，yk,yk+1;h） 它每步求解yk+1需要解一个隐式方程。null欧拉（Euler）方法在x= x0 处，用差商代替导数：由得null同理，在x= xn 处，用差商代替导数：由得若记则上式可记为此即为求解初值问题的Euler方法，又称显式Euler方法。null例: 用Euler方法求解常微分方程初值问题并将数值解和该问题的解析解比较。解：Euler方法的具体格式：nullh=0.2;y(1)=0.2;x=0.2:h:3; for n=1:14 xn=x(n);yn=y(n); y(n+1)=yn+h*(yn/xn-2*yn*yn); end x0=0.2:h:3;y0=x0./(1+x0.^2); plot(x0,y0,x,y,x,y,'o')程序实现null xn y(xn) yn yn-y(xn) 0.0 0 0 0 0.2 0.1923 0.2000 0.0077 0.4 0.3448 0.3840 0.0392 0.6 0.4412 0.5170 0.0758 0.8 0.4878 0.5824 0.0946 1.0 0.5000 0.5924 0.0924 1.2 0.4918 0.5705 0.0787 1.4 0.4730 0.5354 0.0624h=0.2, xn=nh,(n=0,1,2…,15), f(x,y)=y/x – 2y2 计算中取f(0,0)=1. 计算结果如下：nullxn y(xn) yn yn-y(xn) 1.6 0.4494 0.4972 0.0478 1.8 0.4245 0.4605 0.0359 2.0 0.4000 0.4268 0.0268 2.2 0.3767 0.3966 0.0199 2.4 0.3550 0.3698 0.0147 2.6 0.3351 0.3459 0.0108 2.8 0.3167 0.3246 0.0079 3.0 0.3000 0.3057 0.0057由表中数据可以看到，微分方程初值问题的数值解和解析解的误差一般在小数点后第二位或第三位小数上，说明Euler方法的精度是比较差的。null二阶Runge-Kutta方法其中 c1, c2, 2, 21 待定。上式的局部截断误差为：null即c1 = 1-a ， 2 = 21 =1/(2a)方程组解不唯一，可令c2=a  0 ，则 满足上述条件的公式都为2阶R-K公式。null例 蛇形曲线的初值问题令f(x,y)=y/x –2y2, 取 f(0,0)=1, h=0.2,xn=hn , ( n = 1,2,…,15)2阶龙格-库塔公式计算格式： k1=yn/xn – 2yn 2, k2 = (yn+hk1)/(xn+h) – 2(yn+hk1)2 yn+1=yn + 0.5h[ k1 + k2]nullx0=0;y0=0;h=.2; x=.2:h:3; k1=1; k2=(y0+h*k1)/x(1)-2*(y0+h*k1)^2; y(1)=y0+.5*h*(k1+k2);for n=1:14 k1=y(n)/x(n)-2*y(n)^2; k2=(y(n)+h*k1)/x(n+1)-2*(y(n)+h*k1)^2; y(n+1)=y(n)+0.5*h*(k1+k2); endy1=x./(1+x.^2); plot(x,y,'o',x,y1)null常用的一个公式为四阶Runge-Kutta方法null function ydot = harmonic(t,y) ydot=[y(2); -y(1)] y=inline(‘[0 1;-1 0]*y’,’t’,’y’);System of Equationsnull function ydot =twobody(t,y) r=sqrt(y(1)^2+y(2)^2); ydot=[y(3);y(4);-y(1)/r^3;-y(2)/r^3]; Two Body ProblemnullLinearized Differential EquationsnullJ的特征值是解增长解衰减解振荡MATLAB求常微分方程数值解的函数MATLAB求常微分方程数值解的函数 基于龙格－库塔法， MATLAB求常微分方程数值解的函数，一般调用格式为： [t,y]=ode23('fname',tspan,y0) [t,y]=ode45('fname',tspan,y0) 其中fname是定义f(t,y)的函数文件名，该函数文件必须返回一个列向量。tspan形式为[t0,tf],表示求解区间。y0是初始状态列向量。t和y分别给出时间向量和相应的状态向量。null ode23: Bogacki,Shampine(1989)和Shampine(1994), ”23”表示用两算法：一个2阶，一个3阶Bogacki, P. and LF Shampine, "A 3(2) pair of Runge-Kutta formulas," Appl. Math. Letters, Vol. 2, 1989, pp 1-9. BS23 algorithmnull F=inline('[y(2);-y(1)]','t','y') ode23(F,[0 2*pi],[1;0]) opts=odeset(‘reltol’,1.e-4,’abstol’,1.e-6,’outputfcn’)Examplesnull ode23(@twobody,[0 2*pi],[1;0;0;1]);Examplesnull y0=[1;0;0;3]; ode23(@twobody,[0 2*pi],y0); null y0=[1;0;0;3]; [t,y]=ode23(@twobody,[0 2*pi],y0); plot(y(:,1),y(:,2)); axis equalnull y0=[1;0;0;3]; [t,y]=ode23(@twobody,[0 2*pi],y0); plot(y(:,1)) plot(y(:,2))null A problem is stiff if the solution being sought is varying slowly, but there are nearby solutions that very rapidly, so the numerical method must take small steps to obtain satisfactory results.A model of flame propagation（火焰燃烧）: y是球的半径，y^2和y^3与球的表面积和体积有关想一下，点燃一根火柴，光球迅速增长，到达一个关键的大小，然后维持它的大小（由于进入球内氧气和消耗的氧气平衡）Stiff Problem(刚性问题)nulleta=0.02;sym y; F=inline(‘y^2-y^3’,’t’,’y’);ode23(F,[0 2/eta],eta);nulleta=0.00002;sym y; F=inline(‘y^2-y^3’,’t’,’y’);ode23(F,[0 2/eta],eta);nulleta=0.00002; ode23s(inline('y^2-y^3','t','y'),[0 2/eta],eta);null例 蛇形曲线的常微分方程初值问题 MATLAB数值求解命令 F=inline('1./(1+x.^2)-2*y.^2')； ode23(F,[0,6],0) 输出结果为图形 null[T,y]=ode23(f,[0,6],0)将得到自变量和函数的离散数据 T =0 0.0001 0.0005 0.0025 0.0125 0.0496 0.1085 0.1863 0.2837 0.4091 0.5991 0.8513 1.0567 1.2680 1.5110 1.8050 2.1788 2.6842 3.2842 3.8842 4.4842 5.0842 5.6842 6.0000 y =0 0.0001 0.0005 0.0025 0.0125 0.0495 0.1073 0.1800 0.2626 0.3505 0.4411 0.4944 0.5000 0.4868 0.4607 0.4242 0.3793 0.3270 0.2783 0.2411 0.2122 0.1891 0.1705 0.1620null例 洛伦兹模型由如下常微分方程组描述 取 =8/3，=10，=28。 初值：x(0)=0，y(0)=0，z(0)=0.01。 利用MATLAB求解常微分方程数值解命令计算出t∈[0，80]内，三个未知函数的数据值，并绘出相空间在Y-X平面的投影曲线 null气象学家Lorenz提出一篇论文，名叫「一只蝴蝶拍一下翅膀会不会在Taxas州引起龙卷风？」论述某系统如果初期条件差一点点，结果会很不稳定，他把这种现象戏称做「蝴蝶效应」。Lorenz为何要写这篇论文呢？ 这故事发生在1961年的某个冬天，他如往常一般在办公室操作气象电脑。平时，他只需要将温度、湿度、压力等气象数据输入，电脑就会依据三个内建的微分方程式，计算出下一刻可能的气象数据，因此模拟出气象变化图。 这一天，Lorenz想更进一步了解某段纪录的後续变化，他把某时刻的气象数据重新输入电脑，让电脑计算出更多的後续结果。当时，电脑处理数据 资料 新概念英语资料下载李居明饿命改运学pdf成本会计期末资料社会工作导论资料工程结算所需资料清单 的数度不快，在结果出来之前，足够他喝杯咖啡并和友人闲聊一阵。在一小时後，结果出来了，不过令他目瞪口呆。结果和原资讯两相比较，初期数据还差不多，越到後期，数据差异就越大了，就像是不同的两笔资讯。而问题并不出在电脑，问题是他输入的数据差了0.000127，而这些微的差异却造成天壤之别。所以长期的准确预测天气是不可能的。null天气预报的准确性：http://www.nani.com.tw/senior/seniorteach/stfastnews/stnfastnews/stnfastnews4/stnfastnews4_2.htm Lorenz现象的 数学 数学高考答题卡模板高考数学答题卡模板三年级数学混合运算测试卷数学作业设计案例新人教版八年级上数学教学计划 ： http://www.sciscape.org/news_detail.php?news_id=225 分形艺术电子版： http://www.phil.pku.edu.cn/personal/huajie/fractalart/html/notice.htm 混沌理论： http://www.lamost.org/~yzhao/doc/chaos.htmlnull记向量 [y1，y2，y3] = [x，y，z]，创建MATLAB函数文件如下 function z=flo(t,y) z(1,:)=-8*y(1)/3+y(2).*y(3); z(2,:)=-10*(y(2)-y(3)); z(3,:)=-y(1).*y(2)+28*y(2)-y(3);用MATLAB命令求解并绘出Y-X平面的投影图 y0=[0;0;0.01]; [x,y]=ode23(@flo,[0, 80],y0) plot(y(:,2),y(:,1)) nullplot(y(:,3),y(:,1)) nullplot3(y(:,1),y(:,2),y(:,3))null y0=[30;0;-40]; plot(y(:,i))nullnullnull非刚性系统： ode45(Runge-Kutta45) ode23(Runge-Kuatta23) ode113(Adams-Bashforth-Moulton PECE)多步方法刚性系统： ode15s(Gear方法) ,多步方法 ode23s(二阶modified Rosenbroack formula)，单步 ode23t(trapezoidal rule),solve DAEs ode23tb(TR-BDF2) low order methodMatlab's ODE SolversnullLaplacian 算子： Poisson方程（elliptic）：Laplacian 算子的特征值问题:Heat equation(parabolic):Wave equation(hyperbolic):PDE Modelnull五点离散 Poisson方程离散： 特征值问题:Finete Difference Methodsnull热方程：波动方程：nullMatrixsnull椭圆方程：特征值方程：热方程：波动方程：离散方程null 例 用古典显式格式求解抛物型方程 初始条件为 边界条件为 取步长⊿x = h = 0.2 , ⊿t = k = 0.02 。null 解 r = k / h2 = 0.02 / 0.22 = 0.5, 古典显式格式为0 , 0 0.2 , 0 0.4 , 0 1,00 , 0.2 ┇ 0,0.021 , 0.2 ┇ 1, 0.02nullfunction u=gu_dian(f,a,b,c1,c2,m,n) %输入初值和U h=a/(m-1); k=b/(n-1); r=k/h^2; U=zeros(n,m); %赋边界条件 U(2:n,1)=c1; U(2:n,m)=c2;function y=fg(x) y=4.*(x-x.^2);null%赋初始条件 U(1,1:m)=fg(0:h:h*(m-1)); %计算内点上u的数值解U for i=2:n for j=2:(m-1) U(i,j)=(1-2*r)*U(i-1,j)+r*(U(i-1,j-1)+U(i-1,j+1)); end endnull% gu_dianl1.m 步长h=0.20 , k=0.02 , r = k / h2 = 0.5 a=1; b=0.20; c1=0; c2=0; m=6; n=11; U=gu_dian('fg', a,b,c1,c2,m,n) x=0:0.2:a; y=0:0.02:b; [X,Y]=meshgrid(x,y); surf(X,Y,U) % 输入U后再画图nullnullnullPDETOOL有限元方法