3-常微分方程及其求解实例.ppt

第4章 常微分方程及其求解 例： 自变函数 t_final=100; x0=[0;0;1e-10]; % t_final为设定的仿真终止时间 [t,x]=ode45(lorenzeq,[0,t_final],x0); plot(t,x), figure; % 打开新图形窗口 plot3(x(:,1),x(:,2),x(:,3)); axis([10 42 -20 20 -20 25]); % 根据实际数值手动设置坐标系 可采用comet3( )函数绘制动画式的轨迹。 comet3(x(:,1),x(:,2),x(:,3)) 描述微分方程是常微分方程初值问题数值求解的关键。 f1=inline([[-8/3*x(1)+x(2)*x(3); -10*x(2)+10*x(3);... -x(1)*x(2)+28*x(2)-x(3)]],t,x); t_final=100; x0=[0;0;1e-10]; [t,x]=ode45(f1,[0,t_final],x0); plot(t,x), figure; plot3(x(:,1),x(:,2),x(:,3)); axis([10 42 -20 20 -20 25]); 得出完全一致的结果。 （3） 延迟微分方程求解 sol:结构体数据，sol.x:时间向量t, sol.y:状态向量。 例： 编写函数： function dx=c7exdde(t,x,z) xlag1=z(:,1); %第一列表示提取 xlag2=z(:,2); dx=[1-3*x(1)-xlag1(2)-0.2*xlag2(1)^3-xlag2(1); x(3); 4*x(1)-2*x(2)-3*x(3)]; 历史数据函数： function S=c7exhist(t) S=zeros(3,1); function dy=c7exstf2(t,y) dy=[0.04*(1-y(1))-(1-y(2))*y(1)+0.0001*(1-y(2))^2; -10^4*y(1)+3000*(1-y(2))^2]; tic [t2,y2]=ode45(c7exstf2,[0,100],[0;1]); toc length(t2), plot(t2,y2) （5） 隐式微分方程求解 function dx=c7ximp(t,y) A=[sin(y(1)) cos(y(2)); -cos(y(2)) sin(y(1))]; B=[1-y(1); -y(2)]; dx=inv(A)*B; opt=odeset; opt.RelTol=1e-6; [t,y]=ode45(c7ximp,[0,10],[0; 0],opt); plot(t,y) function [x,y]=fdiff(funcs,tspan,x0f,n) t0=tspan(1);tfinal=tspan(2); ga=x0f(1); gb=x0f(2); h=(tfinal-t0)/n; for i=1:n, x(i)=t0+h*(i-1); end, x0=x(1:n-1); t=-2+h^2*feval(funcs,x0,2); tmp=feval(funcs,x0,1); v=1+h*tmp/2; w=1-h*tmp/2; b=h^2*feval(funcs,x0,3); b(1)=b(1)-w(1)*ga; b(n-1)=b(n-1)-v(n-1)*gb; b=b; A=diag(t); for i=1:n-2, A(i,i+1)=v(i); A(i+1,i)=w(i+1); end y=inv(A)*b; x=[x tfinal]; y=[ga; y; gb]; tic a = 0;b = 1; solinit = bvpinit(linspace(a,b,10),[0 0]); sol = bvp4c(@ODEfun,@BCfun,solinit); x = [0:0.02:1];y2 = deval(sol,x); t0=toc tic [t0,y1]=shooting(c7fun1,c7fun2,[0,1],[1;4]); t1=toc tic [t1,y2]=fdiff(c7fun5,[0,1],[1,4],50); t2=toc plot(x,y2(1,:),*,t0,y1(:,1),o,t1,y2,r) IVP应用1: 间歇反应器连串-平行复杂反应系统 在间歇反应器中进行液相反应制备产物B,反应网络如 下,反应可在180~260度温度范围内进行,反应