MATLAB 微分代数方程解法Microsoft Word 文档
第六章用MATLAB求解微分方程及微分方程组
![第六章用MATLAB求解微分方程及微分方程组](https://img.taocdn.com/s3/m/0e6ee56e27d3240c8447ef23.png)
运行结果:u = tan(t-c)
例 2 求微分方程的特解.
d 2 y dy 2 4 29 y 0 dx dx y (0) 0, y ' (0) 15
解 输入命令: y=dsolve('D2y+4*Dy+29*y=0','y(0)=0,Dy(0)=15','x') 运行结果为 : y =3e-2xsin(5x)
不同求解器Solver的特点
求解器 Solver ode23s
OD E类 型
刚 性
特点
说明
刚 ode23tb 性
一步法;2阶 当精度较低时, Rosebrock算法; 计算时间比 低精度 ode15s短 当精度较低时, 梯形算法;低精 计算时间比 度 ode15s短
d 2x dx 2 1000(1 x 2 ) x 0 例 4 dt dt x(0) 2; x' (0) 0 解: 令 y1=x,y2=y1’
y Me k1t 即
模型2:快速饮酒后,体液中酒精含量的变化率
dx k2 x k1Mek1t dt x(0) 0
用Matlab求解模型2:
syms x y k1 k2 M t x=dsolve('Dx+k2*x=k1*M*exp(-k1*t)','x(0)=0','t') 运行结果: M*k1/(-k1+k2)*exp(-k2*t+t*(-k1+k2))-exp(-k2*t)*M*k1/(-k1+k2) 即:
(1)
2. 由于弹头始终对准乙舰,故导弹的速度平行于乙舰与导弹头位置的差向量,
Matlab求解微分方程(组)及偏微分方程(组)
![Matlab求解微分方程(组)及偏微分方程(组)](https://img.taocdn.com/s3/m/82fcfae7763231126fdb1127.png)
第四道 Matlab供解微分圆程(组)之阳早格格创做表里介绍:Matlab供解微分圆程(组)下令供解真例:Matlab供解微分圆程(组)真例本质应用问题通过数教修模所归纳得到的圆程,绝大普遍皆是微分圆程,真真能得到代数圆程的机会很少.另一圆里,不妨供解的微分圆程也是格中有限的,特地是下阶圆程战偏偏微分圆程(组).那便央供咱们必须钻研微分圆程(组)的解法:剖析解法战数值解法.一.相关函数、下令及简介1.正在Matlab中,用大写字母D表示导数,Dy表示y关于自变量的一阶导数,D2y表示y关于自变量的二阶导数,依此类推.函数dsolve用去办理常微分圆程(组)的供解问题,调用要领为:X=dsolve(‘eqn1’,’eqn2’,…)函数dsolve用去解标记常微分圆程、圆程组,如果不初初条件,则供出通解,如果有初初条件,则供出特解.注意,系统缺省的自变量为t2.函数dsolve供解的是常微分圆程的透彻解法,也称为常微分圆程的标记解.然而是,有洪量的常微分圆程虽然从表里上道,其解是存留的,然而咱们却无法供出其剖析解,此时,咱们需要觅供圆程的数值解,正在供常微分圆程数值解圆里,MATLAB 具备歉富的函数,咱们将其统称为solver ,其普遍要领为:[T,Y]=solver(odefun,tspan,y0)证明:(1)solver 为下令ode45、ode23、ode113、ode15s 、ode23s 、ode23t 、ode23tb 、ode15i 之一.(2)odefun 是隐现微分圆程'(,)y f t y =正在积分区间tspan 0[,]f t t =上从0t 到f t 用初初条件0y 供解.(3)如果要赢得微分圆程问题正在其余指定时间面012,,,,f t t t t 上的解,则令tspan 012[,,,]f t t t t =(央供是单调的).(4)果为不一种算法不妨灵验的办理所有的ODE 问题,为此,Matlab 提供了多种供解器solver ,对付于分歧的ODE 问题,采与分歧的solver.表1 Matlab 华文本文献读写函数证明:ode23、ode45是极其时常使用的用去供解非刚刚性的尺度形式的一阶微分圆程(组)的初值问题的解的Matlab时常使用步调,其中:ode23采与龙格-库塔2阶算法,用3阶公式做缺面预计去安排步少,具备矮等的粗度.ode45则采与龙格-库塔4阶算法,用5阶公式做缺面预计去安排步少,具备中等的粗度.3.正在matlab下令窗心、步调或者函数中创修局部函数时,可用内联函数inline,inline函数形式相称于编写M 函数文献,然而不需编写M-文献便不妨形貌出某种数教关系.调用inline函数,只可由一个matlab表白式组成,而且只可返回一个变量,不允许[u,v]那种背量形式.果而,所有央供逻辑运算或者乘法运算以供得最后截止的场合,皆不克不迭应用inline函数,inline函数的普遍形式为:FunctionName=inline(‘函数真质’, ‘所有自变量列表’)比圆:(供解F(x)=x^2*cos(a*x)-b ,a,b是标量;x是背量)正在下令窗心输进:Fofx=inline(‘x .^2*cos(a*x)-b’ , ‘x’,’a’,’b’);g= Fofx([pi/3 pi/3.5],4,1)注意:由于使用内联对付象函数inline 不需要其余修坐m 文献,所有使用比较便当,其余正在使用ode45函数的时间,定义函数往往需要编写一个m 文献去单独定义,那样便当于管造文献,那里不妨使用inline 去定义函数.二.真例介绍例1 供解微分圆程2'2x y xy xe -+= 步调:syms x y; y=dsolve(‘Dy+2*x*y=x*exp(-x^2)’,’x ’) 例2 供微分圆程'0x xy y e +-=正在初初条件(1)2y e =下的特解并画出解函数的图形.步调:syms x y; y=dsolve(‘x*Dy+y-exp(1)=0’,’y(1)=2*exp(1)’,’x ’);ezplot(y)例 3 供解微分圆程组530t dx x y e dt dy x y dt ⎧++=⎪⎪⎨⎪--=⎪⎩正在初初条件00|1,|0t t x y ====下的特解并画出解函数的图形.步调:syms x y t[x,y]=dsolve('Dx+5*x+y=exp(t)','Dy-x-3*y=0','x(0)=1','y(0)=0','t')simple(x);simple(y)ezplot(x,y,[0,1.3]);axis auto2.用ode23、ode45等供解非刚刚性尺度形式的一阶微分圆程(组)的初值问题的数值解(近似解)例4 供解微分圆程初值问题2222(0)1dy y x x dx y ⎧=-++⎪⎨⎪=⎩的数值解,供解范畴为区间[0,0.5].步调:fun=inline('-2*y+2*x^2+2*x','x','y');[x,y]=ode23(fun,[0,0.5],1);plot(x,y,'o-')例5 供解微分圆程22'2(1)0,(0)1,(0)0d y dy y y y y dt dt μ--+===的解,并画出解的图形.12,,7dy x y x dt μ===,则编写M-文献function fy=vdp(t,x)fy=[x(2);7*(1-x(1)^2)*x(2)-x(1)];end正在Matlab 下令窗心编写步调y0=[1;0][t,x]=ode45(@vdp,[0,40],y0);或者[t,x]=ode45('vdp',[0,40],y0);y=x(:,1);dy=x(:,2);plot(t,y,t,dy)训练与思索:M-文献vdp.m 改写成inline 函数步调?Euler 合线法供解的基础思维是将微分圆程初值问题 化成一个代数(好分)圆程,主要步调是用好商()()y x h y x h +-代替微商dydx ,于是 记1,(),k k k k x x h y y x +=+=进而1(),k k y y x h +=+于是例6用Euler 合线法供解微分圆程初值问题的数值解(步少h 与),供解范畴为区间[0,2].分解:本问题的好分圆程为步调:>> clear>> f=sym('y+2*x/y^2');>> a=0;>> b=2;>> h=0.4;>> n=(b-a)/h+1;>> x=0;>> y=1;>> szj=[x,y];%数值解>> for i=1:n-1y=y+h*subs(f,{'x','y'},{x,y});%subs ,替换函数x=x+h;szj=[szj;x,y];end>>szj>> plot(szj(:,1),szj(:,2))证明:替换函数subs比圆:输进subs(a+b,a,4) 意义便是把a用4替换掉,返回4+b,也不妨替换多个变量,比圆:subs(cos(a)+sin(b),{a,b},[sym('alpha'),2])分别用字符alpha替换a战2替换b,返回 cos(alpha)+sin(2)特地证明:本问题可进一步利用四阶Runge-Kutta法供解,Euler合线法本质上便是一阶Runge-Kutta法,Runge-Kutta法的迭代公式为相映的Matlab步调为:>> clear>> f=sym('y+2*x/y^2');>> a=0;>> b=2;>> h=0.4;>> n=(b-a)/h+1;>> x=0;>> y=1;>> szj=[x,y];%数值解>> for i=1:n-1l1=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];end>>szj>> plot(szj(:,1),szj(:,2))训练与思索:(1)ode45供解问题并比较好别.(2)利用Matlab 供微分圆程(4)(3)''20y y y -+=的解. (3)供解微分圆程''2',2(1)0,030,(0)1,(0)0y y y y x y y --+=≤≤==的特解.(4)利用Matlab 供微分圆程初值问题2''''00(1)2,|1,|3x x x y xy y y ==+===的解.指示:尽大概多的思量解法三.微分圆程变更为一阶隐式微分圆程组Matlab 微分圆程解算器只可供解尺度形式的一阶隐式微分圆程(组)问题,果此正在使用ODE 解算器之前,咱们需要干的第一步,也是最要害的一步便是借帮状态变量将微分圆程(组)化成Matlab 可交受的尺度形式.天然,如果ODEs 由一个或者多个下阶微分圆程给出,则咱们应先将它变更成一阶隐式常微分圆程组.底下咱们以二个下阶微分圆程组形成的ODEs 为例介绍怎么样将它变更成一个一阶隐式微分圆程组.Step 1 将微分圆程的最下阶变量移到等式左边,其余移到左边,并按阶次从矮到下排列.形式为:Step 2 为每一阶微分式采用状态变量,最下阶除中注意:ODEs 中所有是果变量的最下阶次之战便是需要的状态变量的个数,最下阶的微分式不需要给它状态变量.Step 3 根据采用的状态变量,写出所有状态变量的一阶微分表白式训练与思索:(1)供解微分圆程组 其中2r =1r =*1,μμ=-1/82.45,μ=(0) 1.2,x =(2)供解隐式微分圆程组提示:使用标记预计函数solve供'''',x y ,而后利用供解微分圆程的要领四.偏偏微分圆程解法Matlab 提供了二种要领办理PDE 问题,一是使用pdepe 函数,它不妨供解普遍的PDEs,具备较大的通用性,然而只收援下令形式调用;二是使用PDE 工具箱,不妨供解特殊PDE 问题,PDEtoll 有较大的限造性,比圆只可供解二阶PDE 问题,而且不克不迭办理片微分圆程组,然而是它提供了GUI 界里,从搀纯的编程中解脱出去,共时还不妨通过File —>Save As 直交死成M 代码.1.普遍偏偏微分圆程(组)的供解(1)Matlab 提供的pdepe 函数,不妨直交供解普遍偏偏微分圆程(组),它的调用要领为:sol=pdepe(m,@pdefun,@pdeic,@pdebc,x,t)@pdefun 是PDE 的问题形貌函数,它必须换成尺度形式:那样,PDE 便不妨编写出心函数:[c,f,s]=pdefun(x,t,u,du),m,x,t 对付应于式中相关参数,du 是u 的一阶导数,由给定的输进变量可表示出c,f,s 那三个函数.@pdebc 是PDE 的鸿沟条件形貌函数,它必须化为形式:于是边值条件不妨编写函数形貌为:[pa,qa,pb,qb]=pdebc(x,t,u,du),其中a 表示下鸿沟,b 表示上鸿沟.@pdeic 是PDE 的初值条件,必须化为形式:00(,)u x t u =,故不妨使用函数形貌为:u0=pdeic(x)sol 是一个三维数组,sol(:,:,i)表示i u 的解,换句话道,k u 对付应x(i)战t(j)时的解为sol(i,j,k),通过sol ,咱们不妨使用pdeval 函数直交预计某个面的函数值.(2)真例证明供解偏偏微分其中, 5.7311.46()x x F x e e -=-且谦脚初初条件12(,0)1,(,0)0u x u x ==及鸿沟条件1(0,)0,u t x ∂=∂221(0,)0,(1,)1,(1,)0u u t u t t x ∂===∂解:(1)对付照给出的偏偏微分圆程战pdepe 函数供解的尺度形式,本圆程改写为 可睹1121220.024()10,,,()10.17u F u u x m c f s F u u u x ∂⎡⎤⎢⎥--⎡⎤⎡⎤∂====⎢⎥⎢⎥⎢⎥-∂⎣⎦⎣⎦⎢⎥⎢⎥∂⎣⎦%目标PDE 函数function [c,f,s]=pdefun(x,t,u,du)c=[1;1];f=[0.024*du(1);0.17*du(2)];temp=u(1)-u(2);s=[-1;1].*(exp(5.73*temp)-exp(-11.46*temp))end(2)鸿沟条件改写为:下鸿沟2010.*00f u ⎡⎤⎡⎤⎡⎤+=⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦上鸿沟1110.*000u f -⎡⎤⎡⎤⎡⎤+=⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦ %鸿沟条件函数function [pa,qa,pb,qb]=pdebc(xa,ua,xb,ub,t)pa=[0;ua(2)];qa=[1;0];pb=[ub(1)-1;0];qb=[0;1];end(3)初值条件改写为:1210u u ⎡⎤⎡⎤=⎢⎥⎢⎥⎣⎦⎣⎦%初值条件函数function u0=pdeic(x)u0=[1;0];end(4)编写主调函数clcx=0:0.05:1;t=0:0.05:2;m=0;sol=pdepe(m,@pdefun,@pdeic,@pdebc,x,t);subplot(2,1,1)surf(x,t,sol(:,:,1))subplot(2,1,2)surf(x,t,sol(:,:,2))训练与思索: This example illustrates the straightforward formulation, computation, and plotting of the solution of a single PDE.This equation holds on an interval 01x ≤≤ for times 0t ≥. The PDE satisfies the initial condition (,0)sin u x x π= and boundary conditions(1)PDEtool (GUI )供解偏偏微分圆程的普遍步调正在Matlab 下令窗心输进pdetool ,回车,PDE 工具箱的图形用户界里(GUI)系统便开用了.从定义一个偏偏微分圆程问题到完毕解偏偏微分圆程的定解,所有历程大概不妨分为六个阶段Step 1 “Draw模式”画造仄里有界地区Ω,通过公式把Matlab系统提供的真体模型:矩形、圆、椭圆战多边形,拉拢起去,死成需要的仄里地区.Step 2 “Boundary模式”定义鸿沟,声明分歧鸿沟段的鸿沟条件.Step 3 “PDE模式”定义偏偏微分圆程,决定圆程典型战圆程系数c,a,f,d,根据简直情况,还不妨正在分歧子地区声明分歧系数.Step 4 “Mesh模式”网格化地区Ω,不妨统造自动死成网格的参数,对付死成的网格举止多次细化,使网格分隔更细更合理.Step 5 “Solve模式”解偏偏微分圆程,对付于椭圆型圆程不妨激活并统造非线性自符合解题器去处理非线性圆程;对付于扔物线型圆程战单直型圆程,树坐初初鸿沟条件后不妨供出给定时刻t的解;对付于特性值问题,不妨供出给定区间上的特性值.供解完毕后,不妨返回到Step 4,对付网格进一步细化,举止再次供解.Step 6 “View模式”预计截止的可视化,不妨通过树坐系统提供的对付话框,隐现所供的解的表面图、网格图、等下线图战箭头梯形图.对付于扔物线型战单直线型问题的解还不妨举止径画演示.(2)真例证明用法供解一个正圆形地区上的特性值问题:正圆形地区为:11,1 1.-≤≤-≤≤x x(1)使用PDE工具箱挨开GUI供解圆程(2)加进Draw模式,画造一个矩形,而后单打矩形,正在弹出的对付话框中树坐Left=-1,Bottom=-1,Width=2,Height=2,确认并关关对付话框(3)加进Boundary模式,鸿沟条件采与Dirichlet条件的默认值(4)加进PDE模式,单打工具栏PDE按钮,正在弹出的对付话框中圆程典型采用Eigenmodes,参数树坐c=1,a=-1/2,d=1,确认后关关对付话框(5)单打工具栏的∆按钮,对付正圆形地区举止初初网格剖分,而后再对付网格进一步细化剖分一次(6)面开solve菜单,单打Parameters选项,正在弹出的对付话框中树坐特性值地区为[-20,20](7)单打Plot菜单的Parameters项,正在弹出的对付话框中选中Color、Height(3-D plot)战show mesh项,而后单打Done确认(8)单打工具栏的“=”按钮,开初供解。
实验二MATLAB 求微分方程的解
![实验二MATLAB 求微分方程的解](https://img.taocdn.com/s3/m/249fcf1aff00bed5b9f31d41.png)
实验二 微分方程求解一、问题背景与实验目的实际应用问题通过数学建模所归纳而得到的方程,绝大多数都是微分方程,真正能得到代数方程的机会很少.另一方面,能够求解的微分方程也是十分有限的,特别是高阶方程和偏微分方程(组).这就要求我们必须研究微分方程(组)的解法,既要研究微分方程(组)的解析解法(精确解),更要研究微分方程(组)的数值解法(近似解).对微分方程(组)的解析解法(精确解),Matlab 有专门的函数可以用,本实验将作一定的介绍.本实验将主要研究微分方程(组)的数值解法(近似解),重点介绍 Euler 折线法.二、相关函数(命令)及简介1.dsolve('equ1','equ2',…):Matlab 求微分方程的解析解.equ1、equ2、…为方程(或条件).写方程(或条件)时用 Dy 表示y 关于自变量的一阶导数,用用 D2y 表示 y 关于自变量的二阶导数,依此类推.2.simplify(s):对表达式 s 使用 maple 的化简规则进行化简. 例如: syms xsimplify(sin(x)^2 + cos(x)^2) ans=13.[r,how]=simple(s):由于 Matlab 提供了多种化简规则,simple 命令就是对表达式 s 用各种规则进行化简,然后用 r 返回最简形式,how 返回形成这种形式所用的规则.例如: syms x[r,how]=simple(cos(x)^2-sin(x)^2) r = cos(2*x) how = combine4.[T,Y] = solver(odefun,tspan,y 0) 求微分方程的数值解. 说明:(1) 其中的 solver 为命令 ode45、ode23、ode113、ode15s 、ode23s 、ode23t 、ode23tb 之一.(2) odefun是显式常微分方程:⎪⎩⎪⎨⎧==00)(),(yt y y t f dt dy(3) 在积分区间 tspan =],[0f t t 上,从0t 到f t ,用初始条件0y 求解.(4) 要获得问题在其他指定时间点 ,210,,t t t 上的解,则令 tspan =],,,[,210f t t t t (要求是单调的).(5) 因为没有一种算法可以有效地解决所有的 ODE 问题,为此,Matlab 提供了多种求解器 Solver ,对于不同的ODE 问题,采用不同的Solver .(6) 要特别的是:ode23、ode45 是极其常用的用来求解非刚性的标准形式的一阶常微分方程(组)的初值问题的解的 Matlab 的常用程序,其中:ode23 采用龙格-库塔2 阶算法,用3 阶公式作误差估计来调节步长,具有低等的精度.ode45 则采用龙格-库塔4 阶算法,用5 阶公式作误差估计来调节步长,具有中等的精度.5.ezplot(x,y ,[tmin,tmax]):符号函数的作图命令.x,y 为关于参数t 的符号函数,[tmin,tmax] 为 t 的取值范围.6.inline():建立一个内联函数.格式:inline('expr', 'var1', 'var2',…) ,注意括号里的表达式要加引号.例:Q = dblquad(inline('y*sin(x)'), pi, 2*pi, 0, pi)三、实验内容1. 几个可以直接用 Matlab 求微分方程精确解的例子: 例1:求解微分方程22xxexy dxdy -=+,并加以验证.求解本问题的Matlab 程序为:syms x y %line1 y=dsolve('Dy+2*x*y=x*exp(-x^2)','x') %line2diff(y ,x)+2*x*y-x*exp(-x^2) %line3 simplify(diff(y ,x)+2*x*y-x*exp(-x^2)) %line4 说明:(1) 行line1是用命令定义x,y 为符号变量.这里可以不写,但为确保正确性,建议写上;(2) 行line2是用命令求出的微分方程的解:1/2*exp(-x^2)*x^2+exp(-x^2)*C1(3) 行line3使用所求得的解.这里是将解代入原微分方程,结果应该为0,但这里给出:-x^3*exp(-x^2)-2*x*exp(-x^2)*C1+2*x*(1/2*exp(-x^2)*x^2+exp(-x^2)*C1)(4) 行line4 用 simplify() 函数对上式进行化简,结果为 0, 表明)(x y y =的确是微分方程的解.例2:求微分方程0'=-+x e y xy 在初始条件e y 2)1(=下的特解,并画出解函数的图形.求解本问题的 Matlab 程序为: syms x yy=dsolve('x*Dy+y-exp(x)=0','y(1)=2*exp(1)','x') ezplot(y)微分方程的特解为:y=1/x*exp(x)+1/x* exp (1) (Matlab 格式),即xe e y x+=,解函数的图形如图 1:图1例3:求微分方程组⎪⎪⎩⎪⎪⎨⎧=--=++035y x dt dy e y x dtdx t在初始条件0|,1|00====t t y x 下的特解,并画出解函数的图形.求解本问题的 Matlab 程序为: syms x y t[x,y]=dsolve('Dx+5*x+y=exp(t)','Dy-x-3*y=0','x(0)=1','y(0)=0','t') simple(x); simple(y);ezplot(x,y ,[0,1.3]);axis auto微分方程的特解(式子特别长)以及解函数的图形均略. 2. 用ode23、ode45等求解非刚性的标准形式的一阶常微分方程(组)的初值问题的数值解(近似解).例4:求解微分方程初值问题⎪⎩⎪⎨⎧=++-=1)0(2222y x x y dxdy 的数值解,求解范围为区间[0, 0.5].fun=inline('-2*y+2*x^2+2*x','x','y'); [x,y]=ode23(fun,[0,0.5],1); x'; y';plot(x,y ,'o-') >> x' ans =0.0000 0.0400 0.0900 0.1400 0.1900 0.2400 0.2900 0.3400 0.3900 0.4400 0.4900 0.5000 >> y' ans =1.0000 0.9247 0.8434 0.7754 0.7199 0.6764 0.6440 0.6222 0.6105 0.6084 0.6154 0.6179 图形结果为图 2.图2例 5:求解描述振荡器的经典的 V er der Pol 微分方程.7,0)0(',1)0(,0)1(222====+--μμy y y dtdy y dty d分析:令,,121dtdx x y x ==则.)1(,1221221x x x dtdx x dtdx --==μ先编写函数文件verderpol.m : function xprime = verderpol(t,x) global mu;xprime = [x(2);mu*(1-x(1)^2)*x(2)-x(1)]; 再编写命令文件vdp1.m : global mu; mu = 7; y0=[1;0][t,x] = ode45('verderpol',[0,40],y0); x1=x(:,1);x2=x(:,2); plot(t,x1)图形结果为图3.图33. 用 Euler 折线法求解前面讲到过,能够求解的微分方程也是十分有限的.下面介绍用 Euler 折线法求微分方程的数值解(近似解)的方法.Euler 折线法求解的基本思想是将微分方程初值问题⎪⎩⎪⎨⎧==00)(),,(yx y y x f dx dy化成一个代数方程,即差分方程,主要步骤是用差商hx y h x y )()(-+替代微商dxdy ,于是:⎪⎩⎪⎨⎧==-+)()),(,()()(00x y y x y x f h x y h x y k k k k 记)(,1k k k k x y y h x x =+=+,从而)(1h x y y k k +=+,则有1,,2,1,0).,(,),(1100-=⎪⎩⎪⎨⎧+=+==++n k y x hf y yh x x x y y k k k k k k 例 6:用 Euler 折线法求解微分方程初值问题⎪⎩⎪⎨⎧=+=1)0(,22y y x y dxdy 的数值解(步长h 取0.4),求解范围为区间[0,2].解:本问题的差分方程为1,,2,1,0).2),( ),(,,4.0,1,021100-=⎪⎪⎪⎩⎪⎪⎪⎨⎧+=+=+====++n k y x y y x f y x hf y y h x x h y x k k k k k k (其中: 相应的Matlab 程序见附录 1. 数据结果为:0 1.0000 0.4000 1.4000 0.8000 2.1233 1.2000 3.1145 1.6000 4.4593 2.0000 6.3074图形结果见图4:图4特别说明:本问题可进一步利用四阶 Runge-Kutta 法求解,读者可将两个结果在一个图中显示,并和精确值比较,看看哪个更“精确”?(相应的 Matlab 程序参见附录 2).四、自己动手1. 求微分方程0sin 2')1(2=-+-x xy y x 的通解.2. 求微分方程x e y y y x sin 5'2''=+-的通解.3. 求微分方程组⎪⎪⎩⎪⎪⎨⎧=-+=++00y x dtdy y x dtdx在初始条件0|,1|00====t t y x 下的特解,并画出解函数()y f x =的图形. 4. 分别用 ode23、ode45 求上述第 3 题中的微分方程初值问题的数值解(近似解),求解区间为[0,2]t ∈.利用画图来比较两种求解器之间的差异.5. 用 Euler 折线法求解微分方程初值问题⎪⎩⎪⎨⎧=-=1)0(,12'32y y xy y 的数值解(步长h 取0.1),求解范围为区间[0,2].6. 用四阶 Runge-Kutta 法求解微分方程初值问题⎩⎨⎧=-=1)0(,cos 'y x e y y x 的数值解(步长h 取0.1),求解范围为区间[0,3].四阶 Runge-Kutta 法的迭代公式为(Euler 折线法实为一阶 Runge-Kutta 法):1,,2,1,0),()2,2()2,2(),()22(6,),(342312143211100-=⎪⎪⎪⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎪⎪⎨⎧++=++=++==++++=+==++n k hL y h x f L L h y h x f L L h y h x f L y x f L L L L L hy y h x x x y y k k k k k k k k k k k k 相应的 Matlab 程序参见附录 2.试用该方法求解第5题中的初值问题. 7. 用 ode45 方法求上述第 6 题的常微分方程初值问题的数值解(近似解),从而利用画图来比较两者间的差异.五、附录附录 1:(fulu1.m)clearf=sym('y+2*x/y^2'); a=0; b=2; h=0.4;n=(b-a)/h+1; x=0; y=1;szj=[x,y]; for i=1:n-1y=y+h*subs(f,{'x','y'},{x,y}); x=x+h;szj=[szj;x,y]; end szjplot(szj(:,1),szj(:,2))附录 2:(fulu2.m)clearf=sym('y-exp(x)*cos(x)'); a=0; b=3; h=0.1;n=(b-a)/h+1; x=0; y=1;szj=[x,y];for i=1:n-1l1=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];endszjplot(szj(:,1),szj(:,2))。
微分方程的matlab求解
![微分方程的matlab求解](https://img.taocdn.com/s3/m/87c812dc5022aaea998f0f7d.png)
数学实验
2、龙格-库塔法
龙格-库塔法是利用泰勒展式将y(x+h)在x处展开,并 取其前面若干项来近似y(x+h)而得到公式 y(x+h) y(x) + h j (x, y(x), h) 如果y(xn) yn,则y(xn+1)的近似值为:
yn+1 = yn + h j (xn, yn, h), n = 0, 1,…
精确解 1 1.0048
1.0187 1.0408 1.0703 1.1065 1.1488 1.1966 1.2493 1.3066 1.3679
上一页 下一页
向前欧拉 向后欧拉 1 1 1 1.0091
matlab求解常微分方程.docx
![matlab求解常微分方程.docx](https://img.taocdn.com/s3/m/9ec62e815ef7ba0d4a733bc1.png)
用matlab求解常微分方程在MATLAB中,由函数dsolve()解决常微分方程(组)的求解问题,其具体格式如下:r 二dsolve('eql,eq2,・••字condl,cond2,・.・;V)匕ql,eq2,・・・*为微分方程或微分方程组,,condl,cond2,.・・;是初始条件或边界条件,P是独立变量,默认的独立变量是讥函数dsolve用来解符号常微分方程、方程组,如果没有初始条件,则求出通解,如果有初始条件,则求出特解。
dy _1例1:求解常微分方程莎一石的MATLAB程序为:dsolve(* Dy=l/(x+y) 1r!x1),注意,系统缺省的自变量为t,因此这里要把自变量写明。
其中:Y=lambertw(X)表示函数关系Y*exp(Y)二X。
例2:求解常微分方程E'-y— 0的MATLAB程序为:Y2=dsolve(1y*D2y-Dy A2=01, 1x f)Y2=dsolve(!D2y*y-Dy A2=0 J )我们看到有两个解,其中一个是常数0。
dx 心 ? —+ 5x + y = e dt空_兀_3『= g2f 例3:求常微分方程组〔力 ' 通解的MATLAB 程序为: [X,Y]=dsolve(f Dx+5*x+y=exp(t),Dy-x-3*y=exp(2*t) 1, 111) [X,Y]=dsolve(1Dx+2 *x-Dy=l0 * cos(t),Dx+Dy+2 *y=4 *exp(-2*t) T ,f x(0)=2,y(0)=0f ,f t T) 以上这些都是常微分方程的精确解法,也称为常微分方程的符号解。
但是,我们知 道,有大量的常微分方程虽然从理论上讲,其解是存在的,但我们却无法求出其解析 解,此吋,我们需要寻求方程的数值解,在求常微分方程数值解方面,MATLAB 具有丰 富的函数,我们将其统称为solver,其一般格式为:i°cosr, 7=2 例4:求常微分方程组 y = 0 z 通解的MATLAB 程序为:[T,Y]=solver(odefun,tspan,yO)该函数表示在区间tspan=[tO,tf]±,用初始条件yO求解显式常微分方程卩=",刃。
Matlab软件求解微分方程_170419
![Matlab软件求解微分方程_170419](https://img.taocdn.com/s3/m/aef5390e866fb84ae45c8d9f.png)
matlab函数定义格式总结matlab中函数定义的一些内容:1, 函数定义格式在matlab中应该做成M文件,文件名要和你文件里的function后面的函数名一致在File新建一个M-file在M-file里编辑函数格式为:function [输出实参表]=函数名(输入实参数)注释部分函数体语句return语句(可以有可以没有)如果是文件中的子函数,则可以任意取名,也可以在同一个文件中定义多个子函数例:function [max,min]=mymainfun(x) %主函数n=length(x);max=mysubfun1(x,n);min=mysubfun2(x);function r=mysubfun1(x,n) %子函数1x1=sort(x);function r=mysubfun2(x) %子函数2x1=sort(x);r=x1(1);Matlab自定义函数的五种方法1、函数文件+调用命令文件:需单独定义一个自定义函数的M文件;2、函数文件+子函数:定义一个具有多个自定义函数的M文件;3、Inline:无需M文件,直接定义;4、Syms+subs:无需M文件,直接定义;5、字符串+subs:无需M文件,直接定义.1、函数文件+调用函数文件:定义多个M文件:%调用函数文件:myfile.mclearclcfor t=1:10y=mylfg(t);fprintf(‘M^(1/3)=%6.4f\n’,t,y);end%自定义函数文件: mylfg.mfunction y=mylfg(x) %注意:函数名(mylfg)必须与文件名(mylfg.m)一致Y=x^(1/3);注:这种方法要求自定义函数必须单独写一个M文件,不能与调用的命令文件写在同一个M文件中。
2、函数文件+子函数:定义一个具有多个子函数的M文件%命令文件:funtry2.mfunction []=funtry2()y=lfg2(t)fprintf(‘M^(1/3)=%6.4f\n’);Endfunction y=lfg2(x)Y= x^(1/3);%注:自定义函数文件funtry2.m中可以定义多个子函数function。
微分代数方程求解matlab
![微分代数方程求解matlab](https://img.taocdn.com/s3/m/a7fb655cf4335a8102d276a20029bd64793e625f.png)
微分代数方程求解matlab微分代数方程是一类常见的数学问题,它涉及到微分方程和代数方程的结合。
在实际应用中,我们经常需要求解这类方程,以获得问题的解析解或数值解。
而MATLAB作为一种强大的数值计算软件,可以帮助我们高效地求解微分代数方程。
首先,让我们来了解一下什么是微分代数方程。
微分代数方程是指同时包含了未知函数及其导数和代数函数的方程。
它可以用来描述许多自然现象和工程问题,如电路、机械系统、生物学模型等。
在MATLAB中,我们可以使用符号计算工具箱来求解微分代数方程。
首先,我们需要定义未知函数及其导数,并将它们表示为符号变量。
然后,我们可以使用符号计算工具箱提供的函数来建立微分代数方程,并求解它们。
例如,考虑一个简单的微分代数方程:dy/dx = x + y。
我们可以使用MATLAB来求解这个方程。
首先,我们定义未知函数y和自变量x为符号变量:syms x y然后,我们建立微分代数方程:eqn = diff(y,x) == x + y接下来,我们可以使用dsolve函数来求解这个微分代数方程:sol = dsolve(eqn)最后,我们可以使用ezplot函数来绘制方程的解析解: ezplot(sol)通过这个简单的例子,我们可以看到MATLAB的强大之处。
它不仅可以帮助我们求解微分代数方程,还可以帮助我们可视化方程的解析解。
除了求解微分代数方程的解析解外,MATLAB还提供了许多数值求解方法。
例如,我们可以使用ode45函数来求解常微分方程组。
这个函数使用了龙格-库塔方法来进行数值积分,并返回方程组的数值解。
总之,MATLAB是一个强大的数值计算软件,可以帮助我们高效地求解微分代数方程。
无论是求解微分代数方程的解析解还是数值解,MATLAB都提供了丰富的工具和函数来满足我们的需求。
通过学习和掌握MATLAB的使用,我们可以更好地理解和应用微分代数方程。
MATLAB微分方程几种求解方法及程序.docx
![MATLAB微分方程几种求解方法及程序.docx](https://img.taocdn.com/s3/m/4cc8378b172ded630b1cb675.png)
第五章控制系统仿真§5.2 微分方程求解方法以一个自由振动系统实例为例进行讨论。
如下图 1 所示弹簧 - 阻尼系统,参数如下:M=5 kg, b=1 N.s/m, k=2 N/m, F=1NxbM Fk图 1 弹簧-阻尼系统假设初始条件为:t时,将m拉向右方,忽略小车的摩擦阻力, x(0) 0m x(0) 0m/ s 求系统的响应。
) 用常微分方程的数值求解函数求解包括ode45、ode23、ode113、ode15s、ode23s 等。
wffc1.m myfun1.m一、常微分方程的数值求解函数ode45 求解解:系统方程为m x b x kx F这是一个单变量二阶常微分方程。
将上式写成一个一阶方程组的形式,这是函数ode45 调用规定的格式。
令:x(1) x(位移)x(2) x x(1)(速度)上式可表示成:x(1)x(2)x(2)x(2)x 1 10* x(2) 20* x(1)下面就可以进行程序的编制。
%写出函数文件 myfun1.mfunction xdot=myfun1(t,x)xdot=[x(2);1-10*x(2)-20*x(1)];%主程序 wffc1.m t=[0 30];x0=[0;0];[tt,xx]=ode45(@myfun1,t,x0);plot(tt,yy(:,1),':b',tt,yy(:,2),'-r')legend('位移','速度')title('微分方程的解x(t)')二、方法 2:m x b x kx FX(s)1G(s)5s2s 2F (s)%用传递函数编程求解ksys1.m num=1;den=[5 1 2];%printsys(num,den)%t=0:0.1:10;sys=tf(num,den);figure(1)step(sys)figure(2)impulse(sys)figure(3)t=[0:0.1:10]';ramp=t;lsim(sys,ramp,t);figure(4)tt=size(t);noise=rand(tt,1);lsim(sys,noise,t)figure(5)yy=0.1*t.^2;lsim(num,den,yy,t)w=logspace(-1,1,100)';[m p]=bode(num,den,w);figure(6)subplot(211);semilogx(w,20*log10(m)); grid onsubplot(212);semilogx(w,p)grid on[gm,pm,wpc,wgc]=margin(sys)figure(7)margin(sys)figure(8)nyquist(sys)figure(9)nichols(sys)方法 3:m x b x kx F5 x x 2x1x 0.2 0.2 x 0.4xInt1Int20.2x''1x'1xs su(t)GsScope0.2G2G10.4x_tTo Workspace Clock%主程序 wffc1.mt=[0 30]; x0=[0;0];[tt,yy]=ode45(@myfun1,t,x0);figure(1)plot(tt,yy(:,1),':b',tt,yy(:,2),'-r')hold onplot(tt,0.2-0.2*yy(:,2)-0.4*yy(:,1),'-.k')legend('位移','速度','加速度')title('微分方程的解 ')figure(2)plot(yy(:,1),yy(:,2))title('平面相轨迹 ')%写出函数文件myfun1.m function xdot=myfun1(t,x)。
Matlab微分方程的解法
![Matlab微分方程的解法](https://img.taocdn.com/s3/m/922e6feae009581b6bd9eb6c.png)
-0.5
-0.55
-0.6
-0.65
-0.7
-0.75
-0.8
-0.85
-0.9
-0.95
-1
0
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
1
time t0=0,tt=1
图3 给定新的初始数据,由函数xprim2定义的ODE解的图形
(d) 求解下面方程组并不难:
x x x x ì ' = - 0.1
在下面的初值问题中,有两个未知函数:x1(t)和x2(t),并用以下式子表达其微... 页码,1/11
Matlab关于微分方程的解法
MATLAB使用龙格-库塔-芬尔格(Runge-Kutta-Fehlberg)方法来解ODE问题。在有限点内计算求解。而 这些点的间距有解的本身来决定。当解比较平滑时,区间内使用的点数少一些,在解变化很快时,区间内应使 用较多的点。 为了得到更多的有关何时使用哪种解法和算法的信息,推荐使用helpdesk。所有求解方程通用的语法或句法在 命令集中头两行给出。时间间隔将以向量t=[t0,tt]给出。 命令ode23可以求解(2,3)阶的常微分方程组,函数ode45使用(4,5)阶的龙格-库塔-芬尔格方法。注意,在这种情 况下x’是x的微分不是x的转置。 在命令集中solver将被诸如ode45函数所取代。
0.6
0.55
0.5
0
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
1
time t0=0,tt=1
图1 由函数xprim1定义的ODE解的图形
(b) 解下面的ODE过程是等价的:
ïíìx' = x2
ïîx(0) = 1
Matlab求解微分方程(组)及偏微分方程(组)
![Matlab求解微分方程(组)及偏微分方程(组)](https://img.taocdn.com/s3/m/b99ecf5eaf1ffc4ffe47acd9.png)
第四讲 Matlab 求解微分方程(组)理论介绍:Matlab 求解微分方程(组)命令 求解实例:Matlab 求解微分方程(组)实例实际应用问题通过数学建模所归纳得到的方程,绝大多数都是微分方程,真正能得到代数方程的机会很少.另一方面,能够求解的微分方程也是十分有限的,特别是高阶方程和偏微分方程(组).这就要求我们必须研究微分方程(组)的解法:解析解法和数值解法. 一.相关函数、命令及简介1.在Matlab 中,用大写字母D 表示导数,Dy 表示y 关于自变量的一阶导数,D2y 表示y 关于自变量的二阶导数,依此类推.函数dsolve 用来解决常微分方程(组)的求解问题,调用格式为:X=dsolve(‘eqn1’,’eqn2’,…)函数dsolve 用来解符号常微分方程、方程组,如果没有初始条件,则求出通解,如果有初始条件,则求出特解.注意,系统缺省的自变量为t2.函数dsolve 求解的是常微分方程的精确解法,也称为常微分方程的符号解.但是,有大量的常微分方程虽然从理论上讲,其解是存在的,但我们却无法求出其解析解,此时,我们需要寻求方程的数值解,在求常微分方程数值解方面,MATLAB 具有丰富的函数,我们将其统称为solver ,其一般格式为:[T,Y]=solver(odefun,tspan,y0)说明:(1)solver 为命令ode45、ode23、ode113、ode15s 、ode23s 、ode23t 、ode23tb 、ode15i 之一.(2)odefun 是显示微分方程'(,)y f t y =在积分区间tspan 0[,]f t t =上从0t 到f t 用初始条件0y 求解.(3)如果要获得微分方程问题在其他指定时间点012,,,,f t t t t 上的解,则令tspan 012[,,,]f t t t t =(要求是单调的).(4)因为没有一种算法可以有效的解决所有的ODE 问题,为此,Matlab 提供了多种求解器solver ,对于不同的ODE 问题,采用不同的solver.表1 Matlab中文本文件读写函数说明:ode23、ode45是极其常用的用来求解非刚性的标准形式的一阶微分方程(组)的初值问题的解的Matlab常用程序,其中:ode23采用龙格-库塔2阶算法,用3阶公式作误差估计来调节步长,具有低等的精度.ode45则采用龙格-库塔4阶算法,用5阶公式作误差估计来调节步长,具有中等的精度.3.在matlab命令窗口、程序或函数中创建局部函数时,可用内联函数inline,inline函数形式相当于编写M函数文件,但不需编写M-文件就可以描述出某种数学关系.调用inline函数,只能由一个matlab表达式组成,并且只能返回一个变量,不允许[u,v]这种向量形式.因而,任何要求逻辑运算或乘法运算以求得最终结果的场合,都不能应用inline函数,inline函数的一般形式为:FunctionName=inline(‘函数内容’, ‘所有自变量列表’)例如:(求解F(x)=x^2*cos(a*x)-b ,a,b是标量;x是向量)在命令窗口输入:Fofx=inline(‘x .^2*cos(a*x)-b’ , ‘x’,’a’,’b’); g= Fofx([pi/3 pi/3.5],4,1) 系统输出为:g=-1.5483 -1.7259注意:由于使用内联对象函数inline 不需要另外建立m 文件,所有使用比较方便,另外在使用ode45函数的时候,定义函数往往需要编辑一个m 文件来单独定义,这样不便于管理文件,这里可以使用inline 来定义函数. 二.实例介绍1.几个可以直接用Matlab 求微分方程精确解的实例 例1 求解微分方程2'2x y xy xe -+=程序:syms x y; y=dsolve(‘Dy+2*x*y=x*exp(-x^2)’,’x ’)例 2 求微分方程'0x xy y e +-=在初始条件(1)2y e =下的特解并画出解函数的图形.程序:syms x y; y=dsolve(‘x*Dy+y-exp(1)=0’,’y(1)=2*exp(1)’,’x ’);ezplot(y)例 3 求解微分方程组530tdx x y e dtdy x y dt⎧++=⎪⎪⎨⎪--=⎪⎩在初始条件00|1,|0t t x y ====下的特解并画出解函数的图形.程序:syms x y t[x,y]=dsolve('Dx+5*x+y=exp(t)','Dy-x-3*y=0','x(0)=1','y(0)=0','t') simple(x); simple(y)ezplot(x,y,[0,1.3]);axis auto2.用ode23、ode45等求解非刚性标准形式的一阶微分方程(组)的初值问题的数值解(近似解)例 4 求解微分方程初值问题2222(0)1dy y x xdx y ⎧=-++⎪⎨⎪=⎩的数值解,求解范围为区间[0,0.5].程序:fun=inline('-2*y+2*x^2+2*x','x','y');[x,y]=ode23(fun,[0,0.5],1); plot(x,y,'o-')例 5 求解微分方程22'2(1)0,(0)1,(0)0d y dyy y y y dt dtμ--+===的解,并画出解的图形.分析:这是一个二阶非线性方程,我们可以通过变换,将二阶方程化为一阶方程组求解.令12,,7dyx y x dtμ===,则 121221212,(0)17(1),(0)0dx x x dtdx x x x x dt⎧==⎪⎪⎨⎪=--=⎪⎩ 编写M-文件vdp.m function fy=vdp(t,x)fy=[x(2);7*(1-x(1)^2)*x(2)-x(1)]; end在Matlab 命令窗口编写程序 y0=[1;0][t,x]=ode45(@vdp,[0,40],y0);或[t,x]=ode45('vdp',[0,40],y0); y=x(:,1);dy=x(:,2); plot(t,y,t,dy)练习与思考:M-文件vdp.m 改写成inline 函数程序? 3.用Euler 折线法求解Euler 折线法求解的基本思想是将微分方程初值问题00(,)()dyf x y dxy x y ⎧=⎪⎨⎪=⎩ 化成一个代数(差分)方程,主要步骤是用差商()()y x h y x h +-替代微商dydx,于是00()()(,())()k k k k y x h y x f x y x h y y x +-⎧=⎪⎨⎪=⎩记1,(),k k k k x x h y y x +=+=从而1(),k k y y x h +=+于是0011(),,0,1,2,,1(,).k k k k k k y y x x x h k n y y hf x y ++=⎧⎪=+=-⎨⎪=+⎩例 6 用Euler 折线法求解微分方程初值问题22(0)1dyx y dxy y ⎧=+⎪⎨⎪=⎩的数值解(步长h 取0.4),求解范围为区间[0,2].分析:本问题的差分方程为00110,1,0.4,0,1,2,,1(,).k k k k k k x y h x x h k n y y hf x y ++===⎧⎪=+=-⎨⎪=+⎩程序:>> clear >> f=sym('y+2*x/y^2'); >> a=0; >> b=2; >> h=0.4; >> n=(b-a)/h+1; >> x=0; >> y=1;>> szj=[x,y];%数值解 >> for i=1:n-1y=y+h*subs(f,{'x','y'},{x,y});%subs ,替换函数 x=x+h; szj=[szj;x,y]; end >>szj>> plot(szj(:,1),szj(:,2))说明:替换函数subs 例如:输入subs(a+b,a,4) 意思就是把a 用4替换掉,返回 4+b ,也可以替换多个变量,例如:subs(cos(a)+sin(b),{a,b},[sym('alpha'),2])分别用字符alpha 替换a 和2替换b ,返回 cos(alpha)+sin(2)特别说明:本问题可进一步利用四阶Runge-Kutta 法求解,Euler 折线法实际上就是一阶Runge-Kutta 法,Runge-Kutta 法的迭代公式为001112341213243(),,(22),6(,),0,1,2,,1(,),22(,),22(,).k k k k k k k k k k k k y y x x x h h y y L L L L L f x y k n h h L f x y L h h L f x y L L f x h y hL ++=⎧⎪=+⎪⎪=++++⎪⎪=⎪=-⎨⎪=++⎪⎪⎪=++⎪⎪=++⎩相应的Matlab 程序为:>> clear >> f=sym('y+2*x/y^2'); >> a=0; >> b=2; >> h=0.4; >> n=(b-a)/h+1; >> x=0; >> y=1;>> szj=[x,y];%数值解 >> for i=1:n-1l1=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]; end>>szj>> plot(szj(:,1),szj(:,2))练习与思考:(1)ode45求解问题并比较差异. (2)利用Matlab 求微分方程(4)(3)''20y y y -+=的解.(3)求解微分方程''2',2(1)0,030,(0)1,(0)0y y y y x y y --+=≤≤==的特解. (4)利用Matlab 求微分方程初值问题2''''00(1)2,|1,|3x x x y xy y y ==+===的解. 提醒:尽可能多的考虑解法 三.微分方程转换为一阶显式微分方程组Matlab 微分方程解算器只能求解标准形式的一阶显式微分方程(组)问题,因此在使用ODE 解算器之前,我们需要做的第一步,也是最重要的一步就是借助状态变量将微分方程(组)化成Matlab 可接受的标准形式.当然,如果ODEs 由一个或多个高阶微分方程给出,则我们应先将它变换成一阶显式常微分方程组.下面我们以两个高阶微分方程组构成的ODEs 为例介绍如何将它变换成一个一阶显式微分方程组.Step 1 将微分方程的最高阶变量移到等式左边,其它移到右边,并按阶次从低到高排列.形式为:()'''(1)'''(1)()'''(1)'''(1)(,,,,,,,,,,)(,,,,,,,,,,)m m n n m n x f t x x x x y y y y y g t x x x x y y y y ----⎧=⎨=⎩Step 2 为每一阶微分式选择状态变量,最高阶除外'''(1)123'''(1)123,,,,,,,,,m m n m m m m n x x x x x x x x x y x y x y x y--++++========注意:ODEs 中所有是因变量的最高阶次之和就是需要的状态变量的个数,最高阶的微分式不需要给它状态变量.Step 3 根据选用的状态变量,写出所有状态变量的一阶微分表达式''''122334123''12123,,,,(,,,,,),,(,,,,,)m m n m m m nm n x x x x x x x f t x x x x xx xg t x x x x +++++======练习与思考:(1)求解微分方程组**'''3312*'''3312()()22x x x y x r r y y y x y r r μμμμμμ⎧+-=+--⎪⎪⎨⎪=+--⎪⎩其中2r =1r =*1,μμ=-1/82.45,μ=(0) 1.2,x =(0)0,y ='(0)0,x ='(0) 1.049355751y =-(2)求解隐式微分方程组''''''''''''2235x y x y x y x y xy y ⎧+=⎨++-=⎩ 提示:使用符号计算函数solve 求'''',x y ,然后利用求解微分方程的方法 四.偏微分方程解法Matlab 提供了两种方法解决PDE 问题,一是使用pdepe 函数,它可以求解一般的PDEs,具有较大的通用性,但只支持命令形式调用;二是使用PDE 工具箱,可以求解特殊PDE 问题,PDEtoll 有较大的局限性,比如只能求解二阶PDE 问题,并且不能解决片微分方程组,但是它提供了GUI 界面,从复杂的编程中解脱出来,同时还可以通过File —>Save As 直接生成M 代码.1.一般偏微分方程(组)的求解(1)Matlab 提供的pdepe 函数,可以直接求解一般偏微分方程(组),它的调用格式为:sol=pdepe(m,@pdefun,@pdeic,@pdebc,x,t)@pdefun 是PDE 的问题描述函数,它必须换成标准形式:(,,)[(,,,)](,,,)m m u u u uc x t x x f x t u s x t u x t x x x-∂∂∂∂∂=+∂∂∂∂∂ 这样,PDE 就可以编写入口函数:[c,f,s]=pdefun(x,t,u,du),m,x,t 对应于式中相关参数,du 是u 的一阶导数,由给定的输入变量可表示出c,f,s 这三个函数.@pdebc 是PDE 的边界条件描述函数,它必须化为形式:(,,)(,,).*(,,,)0up x t u q x t u f x t u x∂==∂ 于是边值条件可以编写函数描述为:[pa,qa,pb,qb]=pdebc(x,t,u,du),其中a 表示下边界,b 表示上边界.@pdeic 是PDE 的初值条件,必须化为形式:00(,)u x t u =,故可以使用函数描述为:u0=pdeic(x)sol 是一个三维数组,sol(:,:,i)表示i u 的解,换句话说,k u 对应x(i)和t(j)时的解为sol(i,j,k),通过sol ,我们可以使用pdeval 函数直接计算某个点的函数值.(2)实例说明 求解偏微分2111222221220.024()0.17()u u F u u t xu u F u u tx ⎧∂∂=--⎪⎪∂∂⎨∂∂⎪=+-⎪∂∂⎩ 其中, 5.7311.46()xx F x e e -=-且满足初始条件12(,0)1,(,0)0u x u x ==及边界条件1(0,)0,u t x ∂=∂221(0,)0,(1,)1,(1,)0uu t u t t x∂===∂ 解:(1)对照给出的偏微分方程和pdepe 函数求解的标准形式,原方程改写为111221220.024()1.*()10.17u u F u u x u F u u u t x x ∂⎡⎤⎢⎥--⎡⎤⎡⎤⎡⎤∂∂∂=+⎢⎥⎢⎥⎢⎥⎢⎥-∂∂∂⎣⎦⎣⎦⎣⎦⎢⎥⎢⎥∂⎣⎦可见1121220.024()10,,,()10.17u F u u x m c f s F u u u x ∂⎡⎤⎢⎥--⎡⎤⎡⎤∂====⎢⎥⎢⎥⎢⎥-∂⎣⎦⎣⎦⎢⎥⎢⎥∂⎣⎦ %目标PDE 函数function [c,f,s]=pdefun(x,t,u,du) c=[1;1];f=[0.024*du(1);0.17*du(2)]; temp=u(1)-u(2);s=[-1;1].*(exp(5.73*temp)-exp(-11.46*temp)) end(2)边界条件改写为:下边界2010.*00f u ⎡⎤⎡⎤⎡⎤+=⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦上边界1110.*000u f -⎡⎤⎡⎤⎡⎤+=⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦%边界条件函数function [pa,qa,pb,qb]=pdebc(xa,ua,xb,ub,t) pa=[0;ua(2)]; qa=[1;0]; pb=[ub(1)-1;0]; qb=[0;1]; end(3)初值条件改写为:1210u u ⎡⎤⎡⎤=⎢⎥⎢⎥⎣⎦⎣⎦%初值条件函数 function u0=pdeic(x) u0=[1;0]; end(4)编写主调函数 clc x=0:0.05:1; t=0:0.05:2; m=0;sol=pdepe(m,@pdefun,@pdeic,@pdebc,x,t); subplot(2,1,1) surf(x,t,sol(:,:,1)) subplot(2,1,2) surf(x,t,sol(:,:,2))练习与思考: This example illustrates the straightforward formulation, computation, and plotting of the solution of a single PDE.2()u u t x xπ∂∂∂=∂∂∂ This equation holds on an interval 01x ≤≤ for times 0t ≥. The PDE satisfies the initial condition (,0)sin u x x π= and boundary conditions(0,)0;(1,)0t uu t e t xπ-∂=+=∂ 2.PDEtool 求解偏微分方程(1)PDEtool (GUI )求解偏微分方程的一般步骤在Matlab 命令窗口输入pdetool ,回车,PDE 工具箱的图形用户界面(GUI)系统就启动了.从定义一个偏微分方程问题到完成解偏微分方程的定解,整个过程大致可以分为六个阶段Step 1 “Draw 模式”绘制平面有界区域Ω,通过公式把Matlab 系统提供的实体模型:矩形、圆、椭圆和多边形,组合起来,生成需要的平面区域.Step 2 “Boundary 模式”定义边界,声明不同边界段的边界条件.Step 3 “PDE 模式”定义偏微分方程,确定方程类型和方程系数c,a,f,d ,根据具体情况,还可以在不同子区域声明不同系数.Step 4 “Mesh 模式”网格化区域Ω,可以控制自动生成网格的参数,对生成的网格进行多次细化,使网格分割更细更合理.Step 5 “Solve 模式”解偏微分方程,对于椭圆型方程可以激活并控制非线性自适应解题器来处理非线性方程;对于抛物线型方程和双曲型方程,设置初始边界条件后可以求出给定时刻t 的解;对于特征值问题,可以求出给定区间上的特征值.求解完成后,可以返回到Step 4,对网格进一步细化,进行再次求解.Step 6 “View 模式”计算结果的可视化,可以通过设置系统提供的对话框,显示所求的解的表面图、网格图、等高线图和箭头梯形图.对于抛物线型和双曲线型问题的解还可以进行动画演示.(2)实例说明用法求解一个正方形区域上的特征值问题:12|0u u u u λ∂Ω⎧-∆-=⎪⎨⎪=⎩ 正方形区域为:11,1 1.x x -≤≤-≤≤(1)使用PDE 工具箱打开GUI 求解方程(2)进入Draw 模式,绘制一个矩形,然后双击矩形,在弹出的对话框中设置Left=-1,Bottom=-1,Width=2,Height=2,确认并关闭对话框(3)进入Boundary 模式,边界条件采用Dirichlet 条件的默认值(4)进入PDE 模式,单击工具栏PDE 按钮,在弹出的对话框中方程类型选择Eigenmodes,参数设置c=1,a=-1/2,d=1,确认后关闭对话框(5)单击工具栏的 按钮,对正方形区域进行初始网格剖分,然后再对网格进一步细化剖分一次(6)点开solve菜单,单击Parameters选项,在弹出的对话框中设置特征值区域为[-20,20](7)单击Plot菜单的Parameters项,在弹出的对话框中选中Color、Height(3-D plot)和show mesh项,然后单击Done确认(8)单击工具栏的“=”按钮,开始求解。
第9章 MATLAB 微分方程问题的求解
![第9章 MATLAB 微分方程问题的求解](https://img.taocdn.com/s3/m/90f45024aaea998fcc220e2b.png)
再用yn代替y(xn),便有欧拉法公式(9.1).
隐式欧拉法和二步欧拉法:
隐式欧拉法 向后差商近似导数
y( x1 )
y ( x1 ) y0 h f ( x1 , y ( x1 ))
y( x1 ) y( x0 ) h
x0
x1
过程,这样的算法称为双步法 , 而前面的二种算法都是单步法 。 由于未知数 yi+1 同时出现在等式的两边,不能直接得到,故 称为隐式 欧拉公式,而前者称为显式 欧拉公式。 二步欧拉法(中点欧拉公式) 一般先用显式 y( x 2 ) y( x0 ) 中心差商近似导数 y( x1 ) 计算一个初值, 2h 再迭代求解。 yi 1 yi 1 2h f ( xi , yi ) i 1, ... , n 1
x0 x1 x2
需要2个初值 y0和 y 来启动递推 y i 1 y i h f ( x i 1 , yi 1 ) ( i 0, ... , n 1) 1
用此公式计算 yi+1 时要用到前两步的信息 yi-1 , yi ,故称为二 步欧拉法(公式)。
Hale Waihona Puke 梯形公式— 显、隐式两种算法的平均 为提高精度,我们再从积分近似的角度来看欧拉公式的 改进。
亚当姆斯显式公式
利用k+1 个节点上的被积函数值 f i , f i 1 , ... , f i k 构造 k 阶牛顿 后插多项式 N k ( xi t h) , t [0, 1], 有
xi 1 xi
f ( x, y( x )) dx N k ( xi t h) h dt Rk ( xi t h) h dt
Step 3: 将 yi+1 与 y( xi+1 ) 在 xi 点的泰勒展开作比较
Matlab解微分方程-25页word资料
![Matlab解微分方程-25页word资料](https://img.taocdn.com/s3/m/39a57b6b59eef8c75ebfb33b.png)
第十六章 偏微分方程的数值解法科学研究和工程技术中的许多问题可建立偏微分方程的数学模型。
包含多个自变量的微分方程称为偏微分方程(partial differential equation),简称PDE 。
偏微分方程问题,其求解是十分困难的。
除少数特殊情况外,绝大多数情况均难以求出精确解。
因此,近似解法就显得更为重要。
本章仅介绍求解各类典型偏微分方程定解问题的差分方法。
16.1 几类偏微分方程的定解问题一个偏微分方程的表示通常如下:(,,,,)x x x y y y x y A B C f x y Φ+Φ+Φ=ΦΦΦ (16.1.1)式中,,,A B C 是常数,称为拟线性(quasilinear)数。
通常,存在3种拟线性方程: 双曲型(hyperbolic)方程:240B AC ->; 抛物线型(parabolic)方程:240B AC -=; 椭圆型(ellliptic)方程:240B AC -<。
16.1.2 双曲型方程最简单形式为一阶双曲型方程:0u ua t x∂∂+=∂∂ (16.1.2) 物理中常见的一维振动与波动问题可用二阶波动方程:22222u u a t x ∂∂=∂∂ (16.1.3) 描述,它是双曲型方程的典型形式。
方程的初值问题为:2222200,(,0)()()t u u at x tx u x x u x x t ϕψ=⎧∂∂=>-∞<<+∞⎪∂∂⎪⎪=⎨⎪∂⎪=-∞<<+∞⎪∂⎩ (16.1.4)边界条件一般有三类,最简单的初边值问题为:2222212000,0(,0)(0,)(),(,)()0()t u ua t T x l t x u x lu t g t u l t g t t T ux x t ϕψ=⎧∂∂==<<<<⎪∂∂⎪⎪=≤⎪⎨==≤≤⎪⎪∂=-∞<<+∞⎪∂⎪⎩ (16.1.5)16.1.3 抛物型方程其最简单的形式为一维热传导方程:220(0)u ua a t x∂∂-=>∂∂ (16.1.8) 方程可以有两种不同类型的定解问题:(1) 初值问题:2200,(,0)()u ua t x t xu x x x ϕ⎧∂∂-=>-∞<<+∞⎪∂∂⎨⎪=-∞<<+∞⎩(16.1.6)(2) 初边值问题:221200,0(,0)()0(0,)(),(,)()0u ua t T x l t x u x x x l u t g t u l t g t t Tϕ⎧∂∂-=<<<<⎪∂∂⎪⎪=≤≤⎨⎪==≤≤⎪⎪⎩(16.1.7) 其中()x ϕ,1()g t ,2()g t 为已知函数,且满足连接条件:12(0)(0),()(0)g l g ϕϕ== (16.1.8)边界条件12(0,)(),(,)()u t g t u l t g t ==为第一类边界条件。
第6章 matlab微分方程问题的解法
![第6章 matlab微分方程问题的解法](https://img.taocdn.com/s3/m/b818e40a79563c1ec5da71a6.png)
例:
>> syms t; u=exp(-5*t)*cos(2*t+1)+5; >> uu=5*diff(u,t,2)+4*diff(u,t)+2*u uu = 87*exp(-5*t)*cos(2*t+1)+92*exp(-5*t)*sin(2*t+1)+10 >> syms t y; >> y=dsolve(['D4y+10*D3y+35*D2y+50*Dy+24*y=',... '87*exp(-5*t)*cos(2*t+1)+92*exp(-5*t)*sin(2*t+1)+10'])
• 求解: >> x0=[1.2; 0; 0; -1.04935751]; >> tic, [t,y]=ode45('apolloeq',[0,20],x0); toc elapsed_time = 0.8310 >> length(t), >> plot(y(:,1),y(:,3)) ans = 689
第六章 微分方程问题的解法
• 微分方程的解析解方法 • 常微分方程问题的数值解法
–微分方程问题算法概述 –四阶定步长 Runge-Kutta算法及 MATLAB 实现 –一阶微分方程组的数值解 –微分方程转换
• 特殊微分方程的数值解 • 边值问题的计算机求解 • 偏微分方程的解
6.1 微分方程的解析解方法
6.2.3 一阶微分方程组的数值解
6.2.3.1 四阶五级Runge-Kutta-Felhberg算法
matlab微分方程组的解法
![matlab微分方程组的解法](https://img.taocdn.com/s3/m/f0403d8388eb172ded630b1c59eef8c75fbf95fe.png)
一、引言1.1 MATLAB在微分方程组求解中的应用MATLAB作为一种强大的数学工具,被广泛应用于微分方程组的求解与模拟分析。
1.2 本文的研究目的和意义本文旨在探讨MATLAB在求解微分方程组方面的应用方法,帮助读者更好地理解和运用MATLAB进行微分方程组的解法,从而提高数学建模和工程仿真的效率与精度。
二、微分方程组的基本概念2.1 微分方程组的定义微分方程组是由多个未知函数及其偏导数构成的方程组。
常见的微分方程组可以分为线性微分方程组与非线性微分方程组。
2.2 微分方程组的求解方法求解微分方程组的方法包括解析解法、数值解法和符号解法。
而MATLAB在微分方程数值解法中具有独特的优势。
三、MATLAB在微分方程组求解中的基本操作3.1 MATLAB中微分方程组的表示在MATLAB中,微分方程组可以使用符号表达式或者函数形式表示,便于进行数值求解和仿真分析。
3.2 MATLAB中微分方程组的数值求解利用MATLAB中的ode45、ode23等求解微分方程组的函数,可以快速地求得微分方程组的数值解,并且可以灵活地控制求解的精度和速度。
3.3 MATLAB中微分方程组的图像绘制MATLAB提供了丰富的绘图函数,能够直观地展现微分方程组的数值解,帮助用户更直观地理解微分方程组的解法结果。
四、 MATLAB在微分方程组求解中的应用实例4.1 简单的线性微分方程组求解通过一个简单的线性微分方程组的求解实例,展示MATLAB在微分方程组求解中的基本操作和方法。
4.2 复杂的非线性微分方程组求解通过一个包含非线性项的微分方程组求解实例,展示MATLAB在处理复杂微分方程组时的应用能力。
五、MATLAB在微分方程组求解中的进阶应用5.1 高阶微分方程组的数值求解MATLAB可以利用符号运算工具箱对高阶微分方程组进行符号求解,也可以通过数值求解的方式得到高阶微分方程组的数值解。
5.2 特定约束条件下的微分方程组求解MATLAB可以通过引入特定的约束条件,对微分方程组进行求解,满足实际应用中的各种约束条件。
MATLAB求解微分方程(实验6)微分方程求解-文档资料
![MATLAB求解微分方程(实验6)微分方程求解-文档资料](https://img.taocdn.com/s3/m/b1d7bb0d87c24028915fc38f.png)
dx
输入:y=dsolve ('Dy=1+y^2') y1=dsolve('Dy=1+y^2','y(0)=1','x')
输出:y= tan(t-C1)
(通解)
y1= tan(x+1/4*pi) (特解)
例2 常系数的二阶微分方程
y' '2 y'3y 0, y(0) 1, y' (0) 0
ode23:组合的2/3阶龙格-库塔-芬尔格算法 ode45:运用组合的4/5阶龙格-库塔-芬尔格算法
用于设定误差限(缺省时设定相对误差10-3, 绝对误差10-6), 命令为:options=odeset(’reltol’,rt,’abstol’,at), rt,at:分别为设定的相对误差和绝对误差.
14
注意: 1、在解n个未知函数的方程组时,x0和x均为n维向量,
m-文件中的待解方程组应以x的分量形式写成。
2、使用Matlab软件求数值解时,高阶微分方程必须 等价地变换成一阶微分方程组。
y(n) f (t, y, y, , y(n1) ) y(0), y(0), , y(n1) (0)
选择一组状态变量
x1 y, x2 y, , xn y(n1)
x1 x2 , x2 x3,
xn f (t, x1, x2 , , xn )
15
注意
x1 x2 , x2 x3,
xn f (t, x1, x2 , , xn )
1、建立M文件函数
function xdot = fun(t,x,y) xdot = [x2(t);x3(t);…;f(t, x1(t), x2(t),…xn(t))];
Matlab解微分方程(ODE+PDE)
![Matlab解微分方程(ODE+PDE)](https://img.taocdn.com/s3/m/761d130fba1aa8114431d9ce.png)
常微分方程:1 ODE解算器简介(ode**)2 微分方程转换3 刚性/非刚性问题(Stiff/Nonstiff)4 隐式微分方程(IDE)5 微分代数方程(DAE)6 延迟微分方程(DDE)7 边值问题(BVP) 偏微分方程(PDEs)Matlab解法偏微分方程:1 一般偏微分方程组(PDEs)的命令行求解2 特殊偏微分方程(PDEs)的PDEtool求解3 陆君安《偏微分方程的MATLAB解法先来认识下常微分方程(ODE)初值问题解算器(solver)[T,Y,TE,YE,IE] = odesolver(odefun,tspan,y0,options)sxint = deval(sol,xint)Matlab中提供了以下解算器:输入参数:odefun:微分方程的Matlab语言描述函数,必须是函数句柄或者字符串,必须写成Matlab规范格式(也就是一阶显示微分方程组),这个具体在后面讲解tspan=[t0 tf]或者[t0,t1,…tf]:微分变量的范围,两者都是根据t0和tf的值自动选择步长,只是前者返回所有计算点的微分值,而后者只返回指定的点的微分值,一定要注意对于后者tspan必须严格单调,还有就是两者数据存储时使用的内存不同(明显前者多),其它没有任何本质的区别y0=[y(0),y’(0),y’’(0)…]:微分方程初值,依次输入所有状态变量的初值,什么是状态变量在后面有介绍options:微分优化参数,是一个结构体,使用odeset可以设置具体参数,详细内容查看帮助输出参数:T:时间列向量,也就是ode**计算微分方程的值的点Y:二维数组,第i列表示第i个状态变量的值,行数与T一致在求解ODE时,我们还会用到deval()函数,deval的作用就是通过结构体solution计算t 对应x值,和polyval之类的很相似!参数格式如下:sol:就是上次调用ode**函数得道的结构体解xint:需要计算的点,可以是标量或者向量,但是必须在tspan范围内该函数的好处就是如果我想知道t=t0时的y值,不需要重新使用ode计算,而直接使用上次计算的得道solution就可以[教程] 微分方程转换为一阶显示微分方程组方法好,上面我们把Matlab中的常微分方程(ODE)的解算器讲解的差不多了,下面我们就具体开始介绍如何使用上面的知识吧!现实总是残酷的,要得到就必须先付出,不可能所有的ODE一拿来就可以直接使用,因此,在使用ODE解算器之前,我们需要做的第一步,也是最重要的一步,借助状态变量将微分方程组化成Matlab可接受的标准形式(一阶显示常微分方程)!如果ODEs由一个或多个高阶微分方程给出,则我们应先将其变换成一阶显式常微分方程组!下面我们以两个高阶微分方程构成的ODEs为例介绍如何将之变换成一个一阶显式常微分方程组。
用Matlab软件求常微分方程的解或通解
![用Matlab软件求常微分方程的解或通解](https://img.taocdn.com/s3/m/035a543a1711cc7931b7168a.png)
《高等数学》实验报告实验人员:系(班): 学号: 姓名: 实验地点:电教楼五号机房实验名称:Matlab 高等数学实验实验时间:2014-6-3 16:30--18:30实验名称:用Matlab 软件求常微分方程的解(或通解)实验目的:熟练掌握Matlab 软件求常微分方程的解(或通解)实验内容:(给出实验程序与运行结果)1、求微分方程的特解. 1、⎪⎩⎪⎨⎧===+-10)0(,6)0(034'22y y y dx dy dx y d 程序:>> dsolve('D2y-4*Dy+3*y','y(0)=6,Dy(0)=10','x')ans = 4*exp(x)+2*exp(3*x)吕梁学院《高等数学》实验报告2、⎪⎩⎪⎨⎧===++0)0(,2)0(044'22y y y dx dy dx y d 程序:>>dsolve('4*D2y+4*Dy+y','y(0)=2,Dy(0)=0','x') ans =2*exp(-1/2*x)+exp(-1/2*x)*x 3、⎪⎩⎪⎨⎧===++15)0(',0)0(029422y y y dx dy dx y d 程序:>>dsolve('D2y+4*Dy+29*y=0','y(0)=9,Dy(0)=15','x') ans = 33/5*exp(-2*x)*sin(5*x)+9*exp(-2*x)*cos(5*x)4、⎪⎩⎪⎨⎧===+-3)0(',0)0(013422y y y dx dy dx y d 程序:>>dsolve('D2y-4*dy+13*y=0','y(0)=0','Dy(0)=3','x') ans = 3/13*sin(13^(1/2)*x)*13^(1/2)-4/13*cos(13^(1/2)*x)*dy+4/13*dy 5、⎪⎩⎪⎨⎧-===--5)0(',0)0(04322y y y dx dy dx y d 程序:>>dsolve('D2y-3*Dy-4*y','y(0)=0,Dy(0)=-5','x') ans = exp(-x)-exp(4*x)二、求齐次非线性微分方程的通解1、133222+=--x y dx dy dx y d 程序:>>dsolve('D2y-2*Dy-3*y=3*x+1','x') ans = exp(-x)*C2+exp(3*x)*C1+1/3-x2、x xe y dx dy dx y d 22265=+- 程序:>>dsolve('D2y-5*Dy+6*y=x*exp(2*x)','x') ans = exp(3*x)*C2+exp(2*x)*C1-1/2*x*exp(2*x)*(2+x)3、x x y dx y d cos 422=+ 程序:>>dsolve('D2y+4*y=x*cos(x)','x') ans = sin(2*x)*C2+cos(2*x)*C1+2/9*sin(x)+1/3*x*cos(x)4、x e y dx y d x cos 22+=+ 程序:>>dsolve('D2y+y=exp(x)','x') ans = sin(x)*C2+cos(x)*C1+1/2*exp(x)>>dsolve('D2y+y=cos(x)','x') ans = sin(x)*C2+cos(x)*C1+1/2*cos(x)+1/2*sin(x)*x 则原式=sin(x)*C2+cos(x)*C1+1/2*exp(x)+sin(x)*C2+cos(x)*C1+1/2*cos(x)+1/2*sin(x)*x 5、x y dx dy dx y d 2sin 5222=+- 程序:>>dsolve('D2y-2*Dy+5*y=sin(2*x)','x') ans =exp(x)*sin(2*x)*C2+exp(x)*cos(2*x)*C1+1/17*sin(2*x)+4/17*cos(2*x)3、微分方程实例1、试求的经过点M (0,1)且在此点与直线相切的积分曲线。
Matlab关于微分方程的解法
![Matlab关于微分方程的解法](https://img.taocdn.com/s3/m/c6b4067001f69e31433294b0.png)
Matlab 关于微分方程的解法MATLAB 使用龙格-库塔-芬尔格(Runge-Kutta-Fehlberg )方法来解ODE 问题。
在有限点内计算求解。
而这些点的间距有解的本身来决定。
当解比较平滑时,区间内使用的点数少一些,在解变化很快时,区间内应使用较多的点。
为了得到更多的有关何时使用哪种解法和算法的信息,推荐使用helpdesk 。
所有求解方程通用的语法或句法在命令集中头两行给出。
时间间隔将以向量t=[t0,tt]给出。
命令ode23可以求解(2,3)阶的常微分方程组,函数ode45使用(4,5)阶的龙格-库塔-芬尔格方法。
注意,在这种情况下x ’是x 的微分不是x 的转置。
在命令集中solver 将被诸如ode45函数所取代。
命令集 龙格-库塔-芬尔格方法[time,x]=solver(str,t,x0) 计算ODE 或由字符串str 给定的ODE 的值,部分解已在向量time 中给出。
在向量time中给出部分解,包含的是时间值。
还有部分解在矩阵x 中给出,x 的列向量是每个方程在这些值下的解。
对于标量问题,方程的解将在向量x 中给出。
这些解在时间区间t(1)到t(2)上计算得到。
其初始值是x0即x(t(1)).此方程组有str 指定的M 文件中函数表示出。
这个函数需要两个参数:标量t 和向量x,应该返回向量x ’(即x 的导数)。
因为对标量ODE 来说,x 和x ’都是标量。
在M 文件中输入odefile 可得到更多信息。
同时可以用命令numjac 来计算Jacobi 函数。
[t,x]=solver(str,t,x0,val) 此方程的求解过程同上,结构val 包含用户给solver 的命令。
参见odeset 和表1,可得到更多信息。
Ode45 此方法被推荐为首选方法。
Ode23 这是一个比ode45低阶的方法。
Ode113 用于更高阶或大的标量计算。
Ode23t 用于解决难度适中的问题。
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
微分代数方程(DAE)的Matlab解法
所谓微分代数方程,是指在微分方程中,某些变量满足某些代数方程的约束。
假设微分方程的更一般形式
可以写成
前面所介绍的ODEs数值解法主要针对能够转换为一阶常微分方程组的类型,故DAE就无法使用前面介绍的常微分方程解法直接求解,必须借助DAE的特殊解法。
其实对于我们使用Matlab求解DAE时,却没有太大的改变只需增加一个Mass参数即可。
描述f(t,x)的方
法和普通微分方程完全一致。
注意:ode15i没法设置Mass属性,换句话说除了ode15i外其他ode计算器都可以求解DAEs问题1.当M(t,y)非奇异的时候,我们可以将微分方程等效的转换为y'=inv(M)*f(t,y),此时就是一个普通的ODE(当
然我们可以将它当成DAEs处理),对任意一个给定的初值条件都有唯一的解
2.当m(t,y)奇异时,我们叫它为DAEs(微分代数方程),DAEs问题只有在同时提供状态变量初值y0和状态变量一阶导数初值py0,且满足M(t0,y0)*yp0=f(t0,y0)时才有唯一解,假如不满足上面的方程,DAEs解算器会将提供的y0和py0作为猜测初始值,并重新计算与提供初值最近的封闭初值
3.质量矩阵可是一个常数矩阵(稀疏矩阵),也可以是一个自定义函数的输出。
但是ode23s只能求解Mass
是常数的DAEs
4.对于Mass奇异的DAEs问题,特别是设置MassSingular为yes时,只能ode15s和ode23t解算器
5.对于DAE我们还有几个参数需要介绍
a.Mass:质量矩阵,不说了,这个是DAE的关键,后面看例子就明白了
b.MStateDependence:质量矩阵M(t,y)是否是y的函数,可以选择none|{weak}|strong,none表示M与
y无关,weak和strong都表示与y相关
c.MvPattern:注意这个必须是稀疏矩阵,S(i,j)=1表示M(t,y)的第i行中任意元素都与第j个状态变量yi有
关,否则为0
d.MassSingular:设置Mass矩阵是否奇异,当设置为yes时,只能使用ode15s和ode23t
e.InitialSlope:状态变量的一阶导数初值yp0,和y0具有相同的size,当使用ode15s和ode23t时,该属
性默认为0
下面我们以实例说明,看下面的例子,求解该方程的数值解
【解】
真是万幸,选取状态变量和求状态变量的一阶导数等,微分方程转换工作,题目已经帮我们完成。
可是细心的网友会发现,最后一个方程不是微分方程而是一个代数方程(这就是为什叫DAE的原因),其实
我们可以将它视为对三个状态变量的约束。
(1)用矩阵形式表示出该DAEs
(2)编写Matlab代码
odefun=@(t,x)[-0.2*x(1)+x(2)*x(3)+0.3*x(1)*x(2);
2*x(1)*x(2)-5*x(2)*x(3)-2*x(2)*x(2);
x(1)+x(2)+x(3)-1];%微分方程组
M=[1 0 0;0 1 0;0 0 0];%质量矩阵
options=odeset('mass',M);%对以DAE问题,mass属性必须设置
x0=[0.8;0.1;0.1];%初值
[t,x]=ode15s(odefun,[0 20],x0,options);%这里好像不能使用ode45 figure('numbertitle','off','name','DAE demo—by Matlabsky')
plot(t,x)
legend('x1(t)','x2(t)','x3(t)')。