Matlab微分方程的解法
重要：MATLAB常微分方程(组)数值解法
Matlab常微分方程求解问题分类
边值问题：
初值问题：
• 定解附加条件在自变量 的一端
• 一般形式为： y' f (x, y)
y(a)
y0
• 初值问题的数值解法一 般采用步进法，如 Runge-Kutta法
➢ 在自变量两端均给定附加 条件
y' f (x, y)
➢ 一般形式：y(a)y1, y(b)y2
1.根据常微分方程要求的求解精度与速度要求
求解初值问题：
y
'
y
2x y
y ( 0 ) 1
(0x1)
比较ode45和ode23的求解精度和速度
ode45和ode23的比较-1
function xODE clear all clc
format long
y0 = 1; [x1,y1] = ode45(@f,[0,1],y0); [x2,y2] = ode23(@f,[0,1],y0); plot(x1,y1,'k-',x2,y2,'b--') xlabel('x') ylabel('y')
rD = k(3)*C(2)-k(5)*C(4);
rE = k(4)*C(3)+k(5)*C(4);
% Mass balances dCdt = [rA; rB; rC; rD; rE];
三个串联的CSTR等温反应器（例4-3）
function IsothermCSTRs clear all clc CA0 = 1.8; % kmol/m^3 CA10 = 0.4; % kmol/m^3 CA20 = 0.2; % kmol/m^3 CA30 = 0.1; % kmol/m^3 k = 0.5; % 1/min tau = 2; stoptime = 2.9; % min [t,y] = ode45(@Equations,[0 stoptime],[CA10 CA20 CA30],[],k,CA0,tau); disp(' Results:') disp(' t CA1 CA2 CA3') disp([t,y]) plot(t,y(:,1),'k--',t,y(:,2),'b:',t,y(:,3),'r-') legend('CA_1','CA_2','CA_3') xlabel('Time (min)') ylabel('Concentration') % -----------------------------------------------------------------function dydt = Equations(t,y,k,CA0,tau) CA1 = y(1); CA2 = y(2); CA3 = y(3); dCA1dt = (CA0-CA1)/tau - k*CA1; dCA2dt = (CA1-CA2)/tau - k*CA2; dCA3dt = (CA2-CA3)/tau - k*CA3; dydt = [dCA1dt; dCA2dt; dCA3dt];
用Matlab解微分方程
用Matlab 软件求解微分方程1．解析解（1）一阶微分方程 求21y dxdy +=的通解：dsolve('Dy=1+y^2','x') 求y x dxdy -+=21的通解：dsolve('Dy=1+x^2-y','x') 求⎪⎩⎪⎨⎧=+=1)0(12y y dx dy 的特解：dsolve('Dy=1+y^2',’y(0)=1’,'x')（2）高阶微分方程 求解⎩⎨⎧-='==-+'+''.2)2(,2)2(,0)(222πππy y y n x y x y x 其中，21=n ，命令为： dsolve('x^2*D2y+x*Dy+(x^2-0.5^2)*y=0','y(pi/2)=2,Dy(pi/2)=-2/pi','x') 求042=-+'-'''x y y y 的通解，命令为：dsolve('D3y-2*Dy+y-4*x=0','x')输出为：ans=8+4*x+C1*exp(x)+C2*exp(-1/2*(5^(1/2)+1)*x)+C3*exp(1/2*(5^(1/2)-1)*x)（3）一阶微分方程组求⎩⎨⎧+-='+=').(3)(4)(),(4)(3)(x g x f x g x g x f x f 的通解：[f,g]=dsolve('Df=3*f+4*g','Dg=-4*f+3*g','x') 输出为： f =exp(3*x)*(cos(4*x)*C1+sin(4*x)*C2)g =-exp(3*x)*(sin(4*x)*C1-cos(4*x)*C2)若再加上初始条件1)0(,0)0(==g f ，则求特解：[f,g]=dsolve('Df=3*f+4*g','Dg=-4*f+3*g','f(0)=0,g(0)=1','x')输出为： f =exp(3*x)*sin(4*x)g =exp(3*x)*cos(4*x)2．数值解（1）一阶微分方程⎪⎩⎪⎨⎧=≤≤-=.1)0(,10,2y x y x y dxdy 现以步长h=0.1用“4阶龙格—库塔公式”求数值解： 先建立“函数M —文件”：function f=eqs1(x,y)f=y-2*x/y;再命令： 格式为：[自变量,因变量]=ode45(‘函数文件名’,节点数组,初始值) 命令为： [x,y]=ode45('eqs1',0:0.1:1,1)若还要画图，就继续命令： plot(x,y)（2）一阶微分方程组⎪⎩⎪⎨⎧==+-='≤≤-+='.3.0)0(,2.0)0(,2sin ,10,2cos 21212211y y y y x y x y y x y 只须向量化，即可用前面方法： function f=eqs2(x,y)f=[cos(x)+2*y(1)-y(2);sin(x)-y(1)+2*y(2)];将此函数文件，以文件名eqs2保存后，再下命令：[x,y]=ode45('eqs2',0:0.1:1,[0.2;0.3])（注：输出的y 是矩阵，第i 列为函数i y 的数值解）要画图，继续命令：hold on,plot(x,y(:,1)),plot(x,y(:,2)),hold off（3）高阶微分方程先化成一阶微分方程组，再用前面方法。
使用Matlab进行微分方程求解的方法
使用Matlab进行微分方程求解的方法引言微分方程是数学中非常重要的一部分，广泛应用于物理、经济、工程等领域。
对于大部分微分方程的解析解往往难以求得，而数值解法则成为了一种常用的解决手段。
Matlab作为一种强大的科学计算软件，也提供了丰富的工具和函数用于求解微分方程，本文将介绍一些常见的使用Matlab进行微分方程求解的方法。
一、数值求解方法1. 欧拉方法欧拉方法是最简单的一种数值求解微分方程的方法，它将微分方程的微分项用差分的方式进行近似。
具体的公式为：y(n+1) = y(n) + hf(x(n), y(n))其中，y(n)表示近似解在第n个点的值，h为步长，f(x, y)为微分方程的右端项。
在Matlab中使用欧拉方法进行求解可以使用ode113函数，通过设定不同的步长，可以得到不同精度的数值解。
2. 中点法中点法是较为精确的一种数值求解微分方程的方法，它的计算公式为：k1 = hf(x(n), y(n))k2 = hf(x(n) + h/2, y(n) + k1/2)y(n+1) = y(n) + k2中点法通过计算两个斜率的平均值来得到下一个点的值，相较于欧拉方法，中点法能提供更精确的数值解。
3. 4阶龙格库塔法龙格库塔法是一类高阶数值求解微分方程的方法，其中4阶龙格库塔法是最常用的一种。
它的计算公式为：k1 = hf(x(n), y(n))k2 = hf(x(n) + h/2, y(n) + k1/2)k3 = hf(x(n) + h/2, y(n) + k2/2)k4 = hf(x(n) + h, y(n) + k3)y(n+1) = y(n) + (k1 + 2k2 + 2k3 + k4)/64阶龙格库塔法通过计算多个斜率的加权平均值来得到下一个点的值，相较于欧拉方法和中点法，它的精度更高。
二、Matlab函数和工具除了可以使用以上的数值方法进行微分方程求解之外，Matlab还提供了一些相关的函数和工具，方便用户进行微分方程的建模和求解。
MATLAB求微分方程的解
实验二微分方程求解一、问题背景与实验目的实际应用问题通过数学建模所归纳而取得的方程，绝大多数都是微分方程，真正能取得代数方程的机遇很少．另一方面，能够求解的微分方程也是十分有限的，专门是高阶方程和偏微分方程（组）．这就要求咱们必需研究微分方程（组）的解法，既要研究微分方程（组）的解析解法（精准解），更要研究微分方程（组）的数值解法（近似解）．对微分方程（组）的解析解法(精准解)，Matlab 有专门的函数能够用，本实验将作必然的介绍．本实验将要紧研究微分方程(组)的数值解法（近似解），重点介绍 Euler 折线法．二、相关函数（命令）及简介1．dsolve('equ1','equ2',…)：Matlab 求微分方程的解析解．equ一、equ 二、…为方程（或条件）．写方程（或条件）时用 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 为命令 ode4五、ode23、ode113、ode15s 、ode23s 、ode23t 、ode23tb 之一．(2) odefun 是显式常微分方程：⎪⎩⎪⎨⎧==00)(),(y t 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：求解微分方程22x xe xy dxdy -=+，并加以验证． 求解本问题的Matlab 程序为：syms x y %line1y=dsolve('Dy+2*x*y=x*exp(-x^2)','x') %line2diff(y,x)+2*x*y-x*exp(-x^2) %line3simplify(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 dtdy e y x dt dx 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,]);axis auto微分方程的特解(式子专门长)和解函数的图形均略．2. 用ode23、ode45等求解非刚性的标准形式的一阶常微分方程(组)的初值问题的数值解(近似解)．例4：求解微分方程初值问题⎪⎩⎪⎨⎧=++-=1)0(2222y x x y dx dy 的数值解，求解范围为区间[0, ]．fun=inline('-2*y+2*x^2+2*x','x','y');[x,y]=ode23(fun,[0,],1);x';y';plot(x,y,'o-')>> x'ans =>> y'ans =图形结果为图 2．图2例 5：求解描述振荡器的经典的 Ver der Pol 微分方程.7,0)0(',1)0(,0)1(222====+--μμy y y dt dy y dty d 分析：令,,121dt dx x y x ==则.)1(,1221221x x x dtdx x dt dx --==μ 先编写函数文件： function xprime = verderpol(t,x)global mu;xprime = [x(2);mu*(1-x(1)^2)*x(2)-x(1)];再编写命令文件：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)(),,(y x y y x f dx dy 化成一个代数方程，即差分方程，要紧步骤是用差商h x 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 y h x x x y y k k k k k k 例 6：用 Euler 折线法求解微分方程初值问题⎪⎩⎪⎨⎧=+=1)0(,22y y x y dx dy的数值解(步长h 取，求解范围为区间[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．附录 1：clearf=sym ('y+2*x/y^2');a=0;b=2;h=;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];endszjplot(szj(:,1),szj(:,2))数据结果为：图形结果见图4：图4专门说明：本问题可进一步利用四阶 Runge-Kutta 法求解，读者可将两个结果在一个图中显示，并和精准值比较，看看哪个更“精准”？（相应的 Matlab 程序参见附录 2）．四、自己动手1. 求微分方程0-x-yx的通解．xy+sin2')1(2=2. 求微分方程x2''='+-的通解．ye5yy x sin3. 求微分方程组⎪⎪⎩⎪⎪⎨⎧=-+=++00y x dtdy y x dt dx 在初始条件0|,1|00====t t y x 下的特解，并画出解函数()y f x =的图形．4. 别离用 ode23、ode45 求上述第 3 题中的微分方程初值问题的数值解(近似解)，求解区间为[0,2]t ∈．利用画图来比较两种求解器之间的不同．5. 用 Euler 折线法求解微分方程初值问题⎪⎩⎪⎨⎧=-=1)0(,12'32y y x y y 的数值解(步长h 取，求解范围为区间[0,2]．6. 用四阶 Runge-Kutta 法求解微分方程初值问题⎩⎨⎧=-=1)0(,cos 'y x e y y x 的数值解(步长h 取，求解范围为区间[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 h y y h x x x y y k k k k k k k k k k k k 相应的 Matlab 程序参见附录 2．附录 2：clearf=sym('y-exp(x)*cos(x)');a=0;b=3;h=;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))试用该方式求解第5题中的初值问题．7. 用 ode45 方式求上述第 6 题的常微分方程初值问题的数值解(近似解)，从而利用画图来比较二者间的不同．五、附录附录 1：clearf=sym('y+2*x/y^2');a=0;b=2;h=;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];endszjplot(szj(:,1),szj(:,2))附录 2：clearf=sym('y-exp(x)*cos(x)');a=0;b=3;h=;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))下面是个例子：% 欲解方程组：% dx/dt=y-x+x(t-1);% dy/dt=y-xy+y(t-1);x0=-1; %初试条件y0=1; %初试条件t=0:01:3;yy=zeros(size(t));xx=zeros(size(t));yy(1)=x0;yy(1)=y0;dt=t(2)-t(1);Nk=round(1/dt);for k=1:length(t)-1;if t(k)<1;k1=yy(k)-xx(k);g1=yy(k)-xx(k)*yy(k);k2=yy(k)+*g1*dt-[xx(k)+*dt*k1];g2=yy(k)+*g1*dt-[xx(k)+*dt*k1]*[yy(k)+*g1*dt];k3=yy(k)+*g2*dt-[xx(k)+*dt*k2];g3=yy(k)+*g2*dt-[xx(k)+*dt*k1]*[yy(k)+*g2*dt];k4=yy(k)+g3*dt-[xx(k)+dt*k3];g4=yy(k)+g3*dt-[xx(k)+dt*k3]*[yy(k)+g3*dt];elsek1=yy(k)-xx(k)+xx(k-Nk);g1=yy(k)-xx(k)*yy(k)+yy(k-Nk);k2=yy(k)+*g1*dt-[xx(k)+*dt*k1]+xx(k-Nk)+*dt*k1;g2=yy(k)+*g1*dt-[xx(k)+*dt*k1]*[yy(k)+*g1*dt]+yy(k-Nk)+*dt*g1;k3=yy(k)+*g2*dt-[xx(k)+*dt*k2]+xx(k-Nk)+*dt*k2;g3=yy(k)+*g2*dt-[xx(k)+*dt*k1]*[yy(k)+*g2*dt]+yy(k-Nk)+*dt*g2;k4=yy(k)+g3*dt-[xx(k)+dt*k3]+xx(k-Nk)+dt*k3;g4=yy(k)+g3*dt-[xx(k)+dt*k3]*[yy(k)+g3*dt]+yy(k-Nk)+dt*k3;endxx(k+1)=xx(k)+dt*[k1+2*k2+2*k3+k4]/6;yy(k+1)=yy(k)+dt*[g1+2*g2+2*g3+g4]/6;endsubplot(121);plot(t,xx);title('t~x(t)');subplot(122);plot(t,yy);title('t~y(t)');/ / / / 文件夹:一、% ４阶龙格库塔法求chen's吸引子% 方程表达式% dx/dt = a*(y-x)% dy/dt = (c-a)*x - x*z + c*y% dz/dt = x*y - b*zclc;close all;clear all;%参数值a = 35;b = 3;c = 28;%初始值x_0 = -1;y_0 = 0;z_0 = 1;h = ; % 积分时间步长step1 = 10000; % 前面的迭代点数step2 = 5000; % 后面的迭代点数X = [];Y = [];Z = [];for(i = 1:1:(step1 + step2))%e1x_e1 = -a*x_0 + a*y_0;y_e1 = c*y_0 + (c - a)*x_0 - x_0*z_0;z_e1 = -b*z_0 + x_0*y_0;%e2y_h = y_0 + *h*y_e1;x_e2 = -a*(x_0 + *h*x_e1) + a*y_h;x_h = x_0 + *h*x_e1;z_h = z_0 + *h*z_e1;y_e2 = c*(y_0 + *h*y_e1) + (c - a)*x_h - x_h*z_h;z_e2 = -b*(z_0 + *h*z_e1) + x_h*y_h;%e3y_h = y_0 + *h*y_e2;x_e3 = -a*(x_0 + *h*x_e2) + a*y_h;x_h = x_0 + *h*x_e2;z_h = z_0 + *h*z_e2;y_e3 = c*(y_0 + *h*y_e2) + (c - a)*x_h - x_h*z_h;z_e3 = -b*(z_0 + *h*z_e2) + x_h*y_h;%e4y_h = y_0 + h*y_e3;x_e4 = -a*(x_0 + h*x_e3) + a*y_h;x_h = x_0 + h*x_e3;z_h = z_0 + h*z_e3;y_e4 = c*(y_0 + h*y_e2) + (c - a)*x_h - x_h*z_h;z_e4 = -b*(z_0 + h*z_e2) + x_h*y_h;%叠代x_1 = x_0 + 1/6*h*(x_e1 + 2*x_e2 +2*x_e3 + x_e4);y_1 = y_0 + 1/6*h*(y_e1 + 2*y_e2 +2*y_e3 + y_e4);z_1 = z_0 + 1/6*h*(z_e1 + 2*z_e2 +2*z_e3 + z_e4);X = [X,x_1];Y = [Y,y_1];Z = [Z,z_1];x_0 = x_1;y_0 = y_1;z_0 = z_1;endX = X(step2+1:end);Y = Y(step2+1:end);Z = Z(step2+1:end);figure(1);plot3(Z(1000:end),Y(1000:end),X(1000:end));grid onxlabel('x');ylabel('y');zlabel('z');这个应该是可以借鉴的。
matlab_常微分方程数值解法
dt 2
简朴问题可以求得解析解，多数实际问题靠数值求解 。
第4页
一阶常微分方程(ODE )初值问题 ： ODE ：Ordinary Differential Equation
dy
f
(x,
y)
dx
x0 x xn
y(x0 ) y0
数值解法就是求y(x)在某些分立旳节点 xn 上旳近似值 yn，用以近似y(xn)
x0
y0
x1 f y(x), x dx
x0
x2 f y(x), x dx
x1
y(x1) f y(x1), x1 h
第17页
同样，在[x0,xn+1] ，积分采用矩形近似，得：
y(xn1) y0
f xn1
x0
y(x), x dx
y(xn ) f y(xn ), xn h
yn y(xn )
第5页
2、欧拉近似办法
2.1 简朴欧拉(L.Euler, 1707-1783)办法。
dy
dx
f
(y, x)
y(x0 ) y0
欧拉数值算法就是由初值通过递推求解，递推求解
就是从初值开始，后一种函数值由前一种函数值得到。核 心是构造递推公式。
y0 y1 y2 yn
第6页
i 1,2,...
第36页
没有一种算法可以有效地解决所有旳 ODE 问题，因此 MATLAB 提供了多种ODE函数。
函数 ODE类
特点
阐明
型
ode45
非刚性 单步法；4，5 阶 R-K 措施；合计 大部分场合旳首选措施
截断误差为 (△x)3
ode23
非刚性 单步法；2，3 阶 R-K 措施；合计 使用于精度较低旳情形
matlab差分法解微分方程
matlab差分法解微分方程在MATLAB中，差分法是一种常用的数值方法，用于解决微分方程。
差分法的基本思想是将微分方程中的导数用离散的差分近似表示，然后通过迭代计算得到方程的数值解。
下面我将从多个角度来解释如何使用差分法在MATLAB中解微分方程。
1. 离散化，首先，我们需要将微分方程离散化，将自变量和因变量分成若干个离散的点。
例如，可以选择一个均匀的网格，将自变量的取值离散化为一系列的点。
这样，微分方程中的导数可以用差分近似来表示。
2. 差分近似，使用差分近似来代替微分方程中的导数。
最常见的差分近似方法是中心差分法。
对于一阶导数，可以使用中心差分公式，f'(x) ≈ (f(x+h) f(x-h)) / (2h)，其中h是离散化步长。
对于二阶导数，可以使用中心差分公式，f''(x) ≈ (f(x+h) 2f(x) + f(x-h)) / (h^2)。
根据微分方程的类型和边界条件，选择适当的差分近似方法。
3. 矩阵表示，将差分近似后的微分方程转化为矩阵形式。
通过将微分方程中的各项离散化，可以得到一个线性方程组。
这个方程组可以用矩阵表示，其中未知量是离散化后的因变量。
4. 数值求解，使用MATLAB中的线性代数求解函数，例如backslash运算符（\）或者LU分解等，求解得到线性方程组的数值解。
这个数值解就是微分方程的近似解。
需要注意的是，差分法是一种数值方法，所得到的解是近似解，精确度受离散化步长的影响。
通常情况下，可以通过减小离散化步长来提高数值解的精确度。
此外，对于某些特殊类型的微分方程，可能需要采用更高级的差分方法，如龙格-库塔法（Runge-Kutta method）或有限元方法（Finite Element Method）等。
综上所述，差分法是一种常用的数值方法，可以在MATLAB中用于解决微分方程。
通过离散化、差分近似、矩阵表示和数值求解等步骤，可以得到微分方程的数值解。
用MATLAB求解微分方程
1. 微分方程的解析解
求微分方程（组）的解析解命令:
dsolve(‘方程1’, ‘方程2’,…‘方程n’, ‘初始条件’, ‘自变量’)
结 果：u = tan(t-c)
解 输入命令：dsolve('Du=1+u^2','t')
STEP2
STEP1
解 输入命令: y=dsolve('D2y+4*Dy+29*y=0','y(0)=0,Dy(0)=15','x')
导弹追踪问题
设位于坐标原点的甲舰向位于x轴上点A(1, 0)处的乙舰发射导弹，导弹头始终对准乙舰.如果乙舰以最大的速度v0(是常数)沿平行于y轴的直线行驶，导弹的速度是5v0，求导弹运行的曲线方程.又乙舰行驶多远时，导弹将它击中？ 解法一（解析法）
由(1),(2)消去t整理得模型:
解法二(数值解)
结 果 为：x = (c1-c2+c3+c2e -3t-c3e-3t)e2t y = -c1e-4t+c2e-4t+c2e-3t-c3e-3t+c1-c2+c3)e2t z = (-c1e-4t+c2e-4t+c1-c2+c3)e2t
2、取t0=0，tf=12，输入命令： [T,Y]=ode45('rigid',[0 12],[0 1 1]); plot(T,Y(:,1),'-',T,Y(:,2),'*',T,Y(:,3),'+')
3、结果如图
图中，y1的图形为实线，y2的图形为“*”线，y3的图形为“+”线.
matlab第6章ode
1 L 0
u
i
状态方程以两个状态元件 i和uo作状态变量 21
+S
i
ui
D
-
dT
L
+
C u0 R
-
0
L
di dt
u0
i
C
d u0
dt
u0
R
T
(2)当S=off：时间长度为: （1-d)*T:
di dt
u0
L
du0
dt
i C
u0
RC
22
S=OFF
u
•
i
•
0
0
1 C
1
1. ode23
在MATLAB中，函数ode23采用2-3阶龙格-库塔法 求解微分方程。
[t,y]=ode23(odefun,tspan,y0) [t,y]=ode23(odefun,tspan,y0,options) odefun：定义微分方程的形式y’=f(t,y) tspan=[t0,tfinal]：表示微分方程的积分限从t0（始值） 到tfinal（终值），该积分限也可以是一些离散的点。 y0:初始状态列向量 options：积分参数，包括‘RelTol’（相对误差）和 ‘AbsTol’（绝对误差），可省略。
IR Vi
L C UC
12
（1）分析：根据电路分析，可以得出微分方程
RI(t) VL(t) VC (t) Vi (t)
I(t) C dVC (t) dt
VL (t)
L
dI(t) dt
LC
d
2VC (t) dt 2
dI(t) L dt RI(t) VC (t) Vi (t)
LC
MATLAB微分方程几种求解方法及程序.docx
第五章控制系统仿真§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微分方程的解法
-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作为一种常用的科学计算软件，在微分方程的数值解法领域具有广泛的应用。
微分方程是描述自然现象中变化规律的数学工具，而数值解法则是指使用计算机进行近似求解微分方程的方法。
在Matlab 中，有多种常用的数值解法可以用来求解微分方程，例如欧拉法、改进的欧拉法、四阶龙格-库塔法等。
本文将对这些数值解法进行介绍和比较，以帮助读者更好地理解和应用微分方程求解数值方法。
二、欧拉法欧拉法是微分方程的最简单的数值解法之一，它通过离散化微分方程进行近似求解。
具体而言，对于一阶常微分方程dy/dx=f(x,y)，可以利用欧拉法进行数值解。
欧拉法的基本思想是将自变量x的增量Δx分成n个小区间，然后根据微分方程的数值近似公式y(x+Δx)=y(x)+f(x,y)Δx对每个小区间进行迭代计算。
欧拉法的优点是简单易实现，但由于它是一阶的数值方法，因此对于某些微分方程求解效果可能不够准确。
三、改进的欧拉法改进的欧拉法是对欧拉法的一种改进，它通过在每个小区间内使用平均斜率来提高求解的精度。
具体而言，对于微分方程dy/dx=f(x,y)，改进的欧拉法可以通过以下迭代公式进行数值求解：y(x+Δx)=y(x)+Δx/2[f(x,y)+f(x+Δx,y+Δx*f(x,y))]改进的欧拉法相比于欧拉法具有更高的数值精度，但计算量也相对增加。
四、四阶龙格-库塔法四阶龙格-库塔法是一种常用的数值微分方程求解方法，它通过四次迭代计算来获得微分方程的数值解。
具体而言，对于微分方程dy/dx=f(x,y)，四阶龙格-库塔法可以用以下公式进行数值求解：k1=f(x,y)k2=f(x+Δx/2,y+Δx/2*k1)k3=f(x+Δx/2,y+Δx/2*k2)k4=f(x+Δx,y+Δx*k3)y(x+Δx)=y(x)+Δx/6*(k1+2*k2+2*k3+k4)四阶龙格-库塔法相比于欧拉法和改进的欧拉法具有更高的数值精度和稳定性，但计算量也相对较大。
matlab求解常微分方程的准确解
matlab求解常微分方程的准确解使用Matlab求解常微分方程的准确解一、引言常微分方程是研究自然界现象和工程实际问题中常见的数学工具之一。
求解常微分方程的准确解对于理解问题的本质和性质具有重要意义。
本文将介绍如何使用Matlab来求解常微分方程的准确解，并通过具体的例子进行演示。
二、常微分方程的基本概念常微分方程是指包含未知函数及其导数的方程。
一般形式为：dy/dx = f(x,y)其中，y是未知函数，x是自变量，f(x,y)是已知函数。
常微分方程的解是指能够满足方程的函数y(x)。
三、Matlab的符号计算工具箱Matlab提供了符号计算工具箱，可以对方程进行符号计算。
通过符号计算工具箱，我们可以求解常微分方程的准确解。
四、使用Matlab求解常微分方程的步骤1. 定义未知函数和自变量。
在Matlab中，可以使用符号变量来定义未知函数和自变量。
2. 定义常微分方程。
使用符号变量来定义常微分方程。
3. 求解常微分方程。
使用dsolve函数来求解常微分方程的准确解。
4. 绘制准确解的图像。
使用ezplot函数来绘制准确解的图像。
五、具体例子假设我们要求解一阶线性常微分方程：dy/dx + y = x其中，y是未知函数，x是自变量。
1. 定义未知函数和自变量。
在Matlab中，可以使用符号变量来定义未知函数和自变量。
syms y(x)2. 定义常微分方程。
使用符号变量来定义常微分方程。
eqn = diff(y,x) + y == x3. 求解常微分方程。
使用dsolve函数来求解常微分方程的准确解。
sol = dsolve(eqn)4. 绘制准确解的图像。
使用ezplot函数来绘制准确解的图像。
ezplot(sol)六、总结本文介绍了如何使用Matlab求解常微分方程的准确解。
通过符号计算工具箱，我们可以方便地求解常微分方程，并得到准确解的图像。
使用Matlab求解常微分方程的准确解可以帮助我们更好地理解问题的本质和性质，并为进一步的分析和应用提供基础。
用Matlab解微分方程
微分方程
程序： t0=0;tf=20; x0=[0,0.25]; [t,x]=ode23(‘f’,t0,tf,x0); plot(t,x(:,1),’:b’,t,x(:,2)，’-r’)
微分方程
• 结果：解函数曲线为如下图形：
微分方程
例3：求y’’=-a2y 的二阶通解: 程序： y=dsolve(‘D2y=-a^2*y’) 结果: y=C1*cos(a*t)+C2*sin(a*t) 例4：求y’’=- a2y y(0)=1,y’(pi/a)=0的二阶特解: 程序： y=dsolve(‘D2y=-a^2*y’,’y(0)=1,Dy(pi/a)=0’) 结果: y=cos(a*t)
微分方程
例5：求微分方程组满足：
df 3 f 4 g dt 初始条件 : f (0) 0, g (0) 1 dg 4 f 3g dt
的解。
微分方程
程序：[f,g]=dsolve(‘Df=3*f+4*g’,’Dg=4*f+3*g’,’f(0)=0,g(0)=1’) 结果： f= exp(3*t)*sin(4*t) g= exp(3*t)*cos(4*t)
用Matlab解微分方程
微分方程
微分方程 一、符号微分方程的求解 dsolve ( ) 用于求解微分方程。 Dy 表示dy/dt (t为缺省的自变量)， Dny 表示y对t的n阶导数。
微分方程
1. 单个方程： 命令格式： y=dsolve(‘Dy=1+y^2’) %求一阶微分方程的通解 y=dsolve(‘Dy=1+y^2’,’y(0)=1’) %求一阶微分方程代初始条件的特解 y=dsolve(‘D2y=cos(2*x),’x’) %求二阶微分方程的通解 y=dsolve(‘D2y=cos(2*x),’y(0)=1’,’Dy(0)=0’,’x’) %求二阶微分方程代初始条件的特解
matlab梯形法解微分方程
主题：matlab梯形法解微分方程内容：一、微分方程的概念和求解方法1. 微分方程的定义2. 微分方程的分类3. 微分方程的解析解和数值解求解方法二、梯形法的原理和步骤1. 梯形法的原理2. 梯形法的求解步骤3. 梯形法的适用范围和优缺点三、matlab中梯形法的实现步骤1. matlab中梯形法的基本函数2. matlab中使用梯形法解微分方程的示例四、实际案例分析1. 利用matlab中的梯形法求解一阶常微分方程2. 利用matlab中的梯形法求解二阶常微分方程五、matlab梯形法解微分方程的应用1. 工程领域中的应用案例2. 科学研究中的应用案例六、总结1. 梯形法解微分方程的优势和局限性2. matlab中梯形法的实际应用效果3. 未来发展方向和展望文章：微分方程是描述自然现象、工程问题等方面中的变化规律的数学工具，它在科学研究和工程应用中都有着重要的地位。
解微分方程的方法有很多种，其中梯形法作为一种数值解方法在matlab中有着丰富的应用。
本文将通过对微分方程的概念、梯形法的原理和步骤、matlab中梯形法的实现步骤、以及实际案例分析，深入探讨matlab梯形法解微分方程的方法和应用。
一、微分方程的概念和求解方法1. 微分方程的定义微分方程是包含一个或多个未知函数及其导数（偏导数）的方程。
根据未知函数、自变量和导数的类型的不同，微分方程可分为常微分方程和偏微分方程。
常微分方程是研究一个未知函数和它的有限阶导数之间的关系的微分方程，而偏微分方程是包含有多个独立变量的方程。
微分方程通常用来描述系统的动力学行为，如弹簧振动、电路的响应等。
2. 微分方程的分类微分方程根据方程中含有未知函数的最高阶导数的阶数、未知函数的个数和自变量的个数等不同特征可以将其分类。
常见的微分方程类型有一阶微分方程、二阶微分方程、线性微分方程、非线性微分方程、常系数微分方程、变系数微分方程等。
3. 微分方程的解析解和数值解求解方法微分方程的解析解法主要包括分离变量法、变参数法、特解法等。
matlab微分方程组的解法
一、引言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解高阶微分方程程序
matlab解高阶微分方程程序（实用版）目录1.引言2.MATLAB 解高阶微分方程的方法3.示例：使用 MATLAB 解一阶微分方程4.示例：使用 MATLAB 解高阶微分方程5.结论正文一、引言微分方程在数学、物理和工程等领域中有着广泛的应用。
在实际问题中，我们常常需要求解高阶微分方程。
使用 MATLAB 这一强大的数学软件工具，可以有效地帮助我们解决这类问题。
本文将介绍如何使用 MATLAB 解高阶微分方程，并给出具体示例。
二、MATLAB 解高阶微分方程的方法在 MATLAB 中，有多种方法可以用来解高阶微分方程。
这里我们将介绍两种常用的方法：符号计算法和数值计算法。
1.符号计算法：MATLAB 的符号计算功能可以用来求解微分方程。
我们可以使用`syms`命令定义符号变量，然后利用`diff`和`int`函数分别求解微分方程的符号解和积分。
2.数值计算法：MATLAB 的数值计算功能可以用来求解微分方程的数值解。
我们可以使用`ode45`、`ode23`等数值积分方法来求解高阶微分方程。
三、示例：使用 MATLAB 解一阶微分方程我们以一阶微分方程为例，演示如何使用 MATLAB 进行符号计算。
假设有一个一阶微分方程：`y" + p(x) * y = q(x)`，其中`p(x)`和`q(x)`是已知函数。
1.使用`syms`命令定义符号变量：`syms(x, y)`。
2.使用`diff`函数求解微分方程的符号解：`dsol = diff(y, x) + p(x) * y = q(x)`。
3.使用`int`函数求解积分：`int(dsol, x)`。
四、示例：使用 MATLAB 解高阶微分方程我们以一个三阶微分方程为例，演示如何使用 MATLAB 进行数值计算。
假设有一个三阶微分方程：`y""" + 2 * y"" + y" + y = f(x)`，其中`f(x)`是已知函数。
二元1阶微分方程的matlab解法
二元1阶微分方程的matlab解法
在MATLAB中，你可以使用dsolve函数来求解二元一阶微分方程。
dsolve函数是MATLAB的符号计算工具箱（Symbolic Math Toolbox）的一部分，它用于求解符号微分方程。
假设你有以下二元一阶微分方程：
dx/dt = f(x, y, t)
dy/dt = g(x, y, t)
其中，f和g是x, y和t的函数。
以下是一个简单的例子，演示如何使用MATLAB的dsolve函数来求解这样的方程。
首先，你需要定义你的微分方程。
在这个例子中，我们假设你有以下方程：
dx/dt = y
dy/dt = -x
这是一个简单的振荡器方程。
然后，你可以使用以下MATLAB代码来求解这个方程：
) t
% 定义微分方程
eq1 = diff(x, t) == y;
eq2 = diff(y, t) == -x;
% 使用dsolve求解
sol = dsolve([eq1, eq2], [x(t), y(t)]);
% 显示解
sol.x
sol.y
导数，然后这些导数被设置为等于f和g。
然后，dsolve函数被用来求解这个方程组，结果存储在sol变量中。
最后，你可以使用sol.x和sol.y来查看x和y的解。
请注意，dsolve函数可能无法求解所有类型的微分方程，特别是那些没有解析解的方程。
对于这些方程，你可能需要使用数值方法（如ode45）来求解。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
有关可用设置
的信息参见表1。
Odeget(set,’set1’)
返回结构set中设置set1的值。
有许多设置对odeset控制的ODE解是有用的，参见表1。例如，如果要在求解过程中画出解的图形，可以 输入：inst=odeset(‘outputfcn’,’odeplot’);.也可使用命令odedemo。
ylabel(‘x values x(0)=-1’); 给出图3
/matlab-pde/ODEmatlab.htm
2004-7-27
在下面的初值问题中，有两个未知函数：x1(t)和x2(t),并用以下式子表达其微... 页码，5/11
x values x(0)=-1
在时间区间t(1)到t(2)上计算得到。其初始值是x0即x(t(1)).此方程组有str指定的
M 文 件 中 函 数 表 示 出。这 个 函 数 需 要 两 个 参 数 ： 标 量 t 和 向 量 x, 应 该 返 回 向 量
x’(即x的导数)。因为对标量ODE来说，x和x’都是标量。在M文件中输入odefile
function xprim=xprim2(t,x) xprim=x.^2; 然后调用MATLAB的ODE算法求解方程。然后画出解的图形： [t,x]=ode45(‘xprim2’,[0,0.95],1); plot(t,x,’o‘,t,x,’-’); xlabel(‘time t0=0,tt=0.95’); ylabel(‘x values x(0)=1’); 得到图2. 注意：在MATLAB中计算出的点在微分绝对值大的区域内更密集些。
表1 ODE求解方程的设置参数
RelTol 给出求解方程允许的相对误差 AbsTol 给出求解方程允许的绝对误差
/matlab-pde/ODEmatlab.htm
2004-7-27
在下面的初值问题中，有两个未知函数：x1(t)和x2(t),并用以下式子表达其微... 页码，2/11
在下面的初值问题中，有两个未知函数：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函数所取代。
/matlab-pde/ODEmatlab.htm
2004-7-27
在下面的初值问题中，有两个未知函数：x1(t)和x2(t),并用以下式子表达其微... 页码，4/11
25
20
x values x(0)=1
15
10
5
0
0
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
'=
1
1 ，这个式子是
以指数形式增长的。大量的被捕食者将会使捕食者的数量增长；同样，越来越少的捕食者会使被捕食者的数量
增长。而且，人口数量也会增长。
创建xprim3,将此函数保存在M文件xprim3.m中：
function xprim=xprim3(t,x) xprim=[x(1)-0.1*x(1)*x(2)+0.01*t;…
用于更高阶或大的标量计算。
Ode23t
用于解决难度适中的问题。
Ode23s Ode15s
用于解决难度较大的微分方程组。对于系统中存在常量矩阵的情况也有用。 与ode23相同，但要求的精度更高。
Ode23tb
用于解决难度较大的问题，对于系统中存在常量矩阵的情况也有用。
Set=odeset(set1,vak1,set2,val2,…) 返回结构set,其中包含用于ODE求解方程的设置参数，
例3 对于某些a和b的值，下面的问题比较难解：
x x x x ì ' = a - (b + 1) + 2
ï1
1
12
x x x x ïïí xï
ï
xïî
' =b
2
0 =1
1
0=3
2
-
1
2 1
2
方程由下面的M文件stiff1.m定义: function stiff=stiff1(t,x) global a; %变量不能放入参数表中 glabol b; stiff=[0,0]; %stiff必须是一个冒号变量 stiff(1)=a-(b+1)*x1+x(1)^2*x(2); stiff(2)=b*x(1)-x(1)^2*x(2); 下面的M文件给出一个比较困难的问题： global a;a=100; global b;b=1; tic; [t,X]=ode23(‘stiff1’，[0 10],[1 3]); toc size(t) 运行后得到的结果如下：
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
首先创建xprim2,将此函数保存在M文件xprim2.m中：
+ 0.01t
ï1
1
12
x x x x ïï
í
' =-
2
+ 0.02
2
1
+ 0.04t
2
ï ïïî
x1 (0) x2 (0)
= =
30 20
x 这个方程组用在人口动力学中。可以认为是单一化的捕食者---被捕食者模式。例如，狐狸和兔子。 1 表示被
x x x 捕食者，
2 表示捕食者。如果被捕食者有无限的食物，并且不会出现捕食者。于是有
120
100
x values x1(0)=30,x2(0)=20
80
60
40
20
x values x1(0)=30,x2(0)=20
0
0
2
4
6
8
10 12 14 16 18 20
time t0=0,tt=20
图4 由函数xprim3定义的ODE解的图形
110
100
90
80
70
60
50
40
30
20
10
-x(2)+0.02*x(1)*x(2)+0.04*t]; 然后调用一个ODE算法和画出解的图形： [t,x]=ode45(‘xprim3’,[0,20],[30;20]);
plot(t,x); xlabel(‘time t0=0,tt=20’); ylabel(‘x values x1(0)=30,x2(0)=20’);
/matlab-pde/ODEmatlab.htm
2004-7-27
在下面的初值问题中，有两个未知函数：x1(t)和x2(t),并用以下式子表达其微... 页码，3/11
x values x(0)=1
1
0.95
0.9
0.85
0.8
0.75
0.7
0.65
可得到更多信息。同时可以用命令numjac来计算Jacobi函数。
[t,x]=solver(str,t,x0,val)
此方程的求解过程同上，结构val包含用户给solver的命令。参见odeset和表1，可得 到更多信息。
Ode45
此方法被推荐为首选方法。
Ode23
这是一个比ode45低阶的方法。
Ode113
Refine 给出与输入点数相乘的因子 OutputFcn 这是一个带有输入函数名的字符串，该字符串将在求解函数执行的每步被调用：odephas2(画出2D
的 平 面 相 位 图 )。Odephas3( 画 出 3D 的 平 面 相 位 图 )，odeplot( 画 出 解 的 图 形 )，odeprint( 显 示 中 间 结 果) OutputSel 是一个整型向量。指出哪些元素应该被传递给函数，特别是传递给OutputFcn Stats 如果参数Stats为on,则将统计并显示出计算过程中资源消耗情况 Jacobian… 如果编写ODE文件代码以便F（t,y,jocobian）返回dF/dy,则将jacobian设置为on Jconstant…如果雅可比数df/dy是常量，则将此参数设置为on Jpattern… 如果编写ODE文件的编码以便函数F（[],[],jpattern）返回带有零的稀疏矩阵并输出非零元素dF/fy,则 需将Jpattern设置为on Vectorized…如果编写ODE文件的编码以便函数F(t,[y1,y2……])返回[F(t,y1) F(t,y2)…],则将此参数设置成on Events… 如果ODE文件中带有参数‘events’，则将此参数设置成on Mass… 如果编写ODE文件编码以实现函数F(t,[],‘mass’)返回M和M(t),应将此参数设置成on MassConstant…如果矩阵M(t)是常量，则将此参数设置成on MaxStep… 此参数是限定算法能使用的区间长度上限的标量 InitialStep… 给出初始步长的标量。如果给定的区间太大，算法就使用一个较小的步长 MaxOrder… 此参数只能被ode15s使用，它主要是指定ode15s的最高阶数,并且此参数应是从1到5的整数 BDF… 此参数只能被ode15s使用，如果倒推微分公式而不是使用通常所使用的微分公式，则要将它设置为