MATLAB解微分方程new.ppt

MATLAB ODE初值问题的数值解 PDE问题的数值解 数值解法需要把连续性的问题加以离散化，从而求出离散节点的数值解。 数值解和精确解 用数值方法求解初值问题，不是求出它的解析解或其近似解析式，而是给出它的解在某些离散节点上的近似值。 单步法和多步法 单步法：在计算yi+1 时只利用y i 多步法：在计算yi+1 时不仅利用y i , 还要利用 yi?1, yi?2,…, xn y(xn) yn yn-y(xn) 0.0 0 0 0 0.2 0.1923 0.2000 0.0077 0.4 0.3448 0.3840 0.0392 0.6 0.4412 0.5170 0.0758 0.8 0.4878 0.5824 0.0946 1.0 0.5000 0.5924 0.0924 1.2 0.4918 0.5705 0.0787 1.4 0.4730 0.5354 0.0624 xn y(xn) yn yn-y(xn) 1.6 0.4494 0.4972 0.0478 1.8 0.4245 0.4605 0.0359 2.0 0.4000 0.4268 0.0268 2.2 0.3767 0.3966 0.0199 2.4 0.3550 0.3698 0.0147 2.6 0.3351 0.3459 0.0108 2.8 0.3167 0.3246 0.0079 3.0 0.3000 0.3057 0.0057 即 MATLAB求常微分方程数值解的函数 基于龙格－库塔法， MATLAB求常微分方程数值解的函数，一般调用格式为： [t,y]=ode23(fname,tspan,y0) [t,y]=ode45(fname,tspan,y0) 其中fname是定义f(t,y)的函数文件名，该函数文件必须返回一个列向量。tspan形式为[t0,tf],表示求解区间。y0是初始状态列向量。t和y分别给出时间向量和相应的状态向量。 例 用古典显式格式求解抛物型方程 初始条件为 边界条件为 取步长⊿x = h = 0.2 , ⊿t = k = 0.02 。 解 r = k / h2 = 0.02 / 0.22 = 0.5, 古典显式格式为 function u=gu_dian(f,a,b,c1,c2,m,n) %输入初值和U h=a/(m-1); k=b/(n-1); r=k/h^2; U=zeros(n,m); %赋边界条件 U(2:n,1)=c1; U(2:n,m)=c2; %赋初始条件 U(1,1:m)=fg(0:h:h*(m-1)); %计算内点上u的数值解U for i=2:n for j=2:(m-1) U(i,j)=(1-2*r)*U(i-1,j)+r*(U(i-1,j-1)+U(i-1,j+1)); end end % gu_dianl1.m 步长h=0.20 , k=0.02 , r = k / h2 = 0.5 a=1; b=0.20; c1=0; c2=0; m=6; n=11; U=gu_dian(fg, a,b,c1,c2,m,n) x=0:0.2:a; y=0:0.02:b; [X,Y]=meshgrid(x,y); surf(X,Y,U) % 输入U后再画图 F=inline([y(2);-y(1)],t,y) ode23(F,[0 2*pi],[1;0]) opts=odeset(‘reltol’,1.e-4,’abstol’,1.e-6,’outputfcn’) Examples ode23(@twobody,[0 2*pi],[1;0;0;1]); 0 1 2 3 4 5 6 7 -1.5 -1 -0.5 0 0.5 1 1.5 Examples y0=[1;0;0;3]; ode23(@twobody,[0 2*pi],y0); 0 1 2 3 4 5 6 7 -2 0 2 4 6 8 10 12 14 16 18 y0=[1;0;0;3]; [t,y]=ode23(@twobody,[0 2*pi],y0); plot(y(:,1),y(:,2)); axis equal y0=[1;0;0;3]; [t,y]=ode23(@twobody,[0 2*pi],y0); plot(y(:,1)) plot(y(:,2)) A problem is stiff if the solution being sought is varying slowly, but there are nearby