11.微分方程问题的解法.ppt

例：椭圆微分方程的拉普拉斯形式： 其自变量取值范围D为： 边界条件为： 1、调整坐标范围:options-Axes Limits 2、设定矩形区域，双击确定调整。 3、设置边界条件：单击边界按钮 ，并双击相应的边弹出对话框，设定边界条件。 4、设定方程。按“PDE”按钮选定Elliptic. 5、单击 按纽可将求解区域用三角形划分成网格。可用 按纽加密。 6、单击三维图标，设置所画曲线特性。 7、单击“plot”或工具栏中的“=”按钮。 例： 编写函数???： 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); 求解： lags=[1 0.5]; tx=dde23(c7exdde,lags,zeros(3,1),[0,10]); plot(tx.x,tx.y(2,:)) ％与ode45()等返回的x矩阵不一样，它是按行排列的。 11.4边值问题的计算机求解 11.4.1 边值问题的打靶算法 数学方法描述： 以二阶方程为例 编写函数： 线性的 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=c7fun1(t,x) xdot=[x(2); -2*x(1)+3*x(2)]; function xdot=c7fun2(t,x) xdot=[x(2); t-2*x(1)+3*x(2)]; [t,y]=shooting(c7fun1, … c7fun2,[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 非线性方程边值问题的打靶算法： 用Newton迭代法处理 m 编写函数： 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(end,3)); end [t,y]=ode45(funcn,tspan,[ga;m],varargin); 例： 编写两个函数： function xdot=c7fun3(t,x) xdot=[x(2); 2*x(1)*x(2); x(4); 2*x(2)*x(3)+2*x(1)*x(4)]; function xdot=c7fun4(t,x) xdot=[x(2); 2*x(1)*x(2)]; [t,y]=nlbound(c7fun4,c7fun3,[0,pi/2],[-1,1],1e-8); plot(t,y); set(gca,xlim,[0,pi/2]); 精确解： 检验： y0=tan(t-pi/4); norm(y(:,1)-y0) ans = 1.66