科学工程计算与matlab编程5-4.pptVIP

5.3.3 微分代数方程求解 编写函数 function dx=c5eqdae(t,x) dx=[-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; options.Mass=M; ％ Mass微分代数方程中 的质量矩阵（控制参数） x0=[0.8; 0.1; 0.1]; [t,x]=ode15s(@c5eqdae,[0,20],x0,options); plot(t,x) 编写函数： function dx=c5eqdae1(t,x) dx=[-0.2*x(1)+x(2)*(1-x(1)-x(2))+0.3*x(1)*x(2); 2*x(1)*x(2)-5*x(2)*(1-x(1)-x(2))-2*x(2)*x(2)]; x0=[0.8; 0.1]; fDae=inline([[-0.2*x(1)+x(2)*(1-x(1)-x(2))+0.3*x(1)*x(2);,... 2*x(1)*x(2)-5*x(2)*(1-x(1)-x(2))-2*x(2)*x(2)]],t,x); [t1,x1]=ode45(fDae,[0,20],x0); plot(t1,x1,t1,1-sum(x1)) 5.3.4延迟微分方程求解 sol:结构体数据，sol.x:时间向量t, sol.y:状态向量。 例： 编写函数： function dx=c5exdde(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)]; 求解： lags=[1 0.5]; tx=dde23(c5exdde,lags,zeros(3,1),[0,10]); plot(tx.x,tx.y(2,:)) ％与ode45()等返回的x矩阵不一样，它是按行排列的。 5.4边值问题的计算机求解 数学方法描述（以二阶方程为例） 编写函数： 线性的 function [t,y]=shooting(f1,f2,tspan,x0f,varargin) t0=tspan(1); tfinal=tspan(2); ga=x0f(1); gb=x0f(2); [t,y1]=ode45(f1,tspan,[1;0],varargin); [t,y2]=ode45(f1,tspan,[0;1],varargin); [t,yp]=ode45(f2,tspan,[0;0],varargin); m=(gb-ga*y1(end,1)-yp(end,1))/y2(end,1); [t,y]=ode45(f2,tspan,[ga;m],varargin); 例： 编写函数： function xdot=c5fun1(t,x) xdot=[x(2); -2*x(1)+3*x(2)]; function xdot=c5fun2(t,x) xdot=[x(2); t-2*x(1)+3*x(2)]; [t,y]=shooting(c5fun1, … c5fun2,[0,1],[1;2]); plot(t,y) 原方程的解析解为 解的检验 y0=((exp(2)-3)*exp(t)+(3-exp(1))*exp(2*t))/(4*exp(1)*(exp(1)-1))+3/4+t/2; norm(y(:,1)-y0) % 整个解函数检验 ans = 4.4790e-008 norm(y(end,1)-2) % 终点条件检验 ans = 2.2620e-008 非线性方程边值问题的打靶算法： 编写函数： function [t,y]=nlbound(funcn,funcv,tspan,x0f,tol,varargin) t0=tspan(1);tfinal=tspan(2); ga=x0f(1); gb=x0f(2); m=1; m0=0; while (norm(m-m0)tol), m0=m; [t,v]=ode45(funcv,tspan,[ga;m;0;1],varargin); m=m0-(v(end,1)-gb)/(v(