资源预览内容
第1页 / 共5页
第2页 / 共5页
第3页 / 共5页
第4页 / 共5页
第5页 / 共5页
亲,该文档总共5页全部预览完了,如果喜欢就下载吧!
资源描述
对流方程偏微分方程的数值解法用迎风格式解对流方程function u = peYF(a,dt,n,minx,maxx,M) format long; h = (maxx-minx)/(n-1); if a0for j=1:(n+M)u0(j) = IniU(minx+(j-M-1)*h);end else for j=1:(n+M)u0(j) = IniU(minx+(j-1)*h);end end u1 = u0; for k=1:Mif a0 for i=(k+1):n+Mu1(i) = -dt*a*(u0(i)-u0(i-1)/h+u0(i);endelsefor i=1:n+M-ku1(i) = -dt*a*(u0(i+1)-u0(i)/h+u0(i);endendu0 = u1; end if a0u = u1(M+1):M+n); elseu = u1(1:n); end format long; 用拉克斯-弗里德里希斯格式解对流方程function u = peHypbLax(a,dt,n,minx,maxx,M) format long; h = (maxx-minx)/(n-1); for j=1:(n+2*M)u0(j) = IniU(minx+(j-M-1)*h); end u1 = u0; for k=1:Mfor i=k+1:n+2*M-ku1(i) = -dt*a*(u0(i+1)-u0(i-1)/h/2+(u0(i+1)+u0(i-1)/2;endu0 = u1; end u = u1(M+1):(M+n); format short; 用拉克斯-温德洛夫格式解对流方程function u = peLaxW(a,dt,n,minx,maxx,M) format long; h = (maxx-minx)/(n-1); for j=1:(n+2*M)u0(j) = IniU(minx+(j-M-1)*h); end u1 = u0; for k=1:Mfor i=k+1:n+2*M-ku1(i) = dt*dt*a*a*(u0(i+1)-2*u0(i)+u0(i-1)/2/h/h - .dt*a*(u0(i+1)-u0(i-1)/h/2+u0(i);endu0 = u1; end u = u1(M+1):(M+n); format short; 用比姆-沃明格式解对流方程function u = peBW(a,dt,n,minx,maxx,M) format long; h = (maxx-minx)/(n-1); for j=1:(n+2*M)u0(j) = IniU(minx+(j-2*M-1)*h); end u1 = u0; for k=1:Mfor i=2*k+1:n+2*Mu1(i) = u0(i)-dt*a*(u0(i)-u0(i-1)/h-a*dt*(1-a*dt/h)* .(u0(i)-2*u0(i-1)+u0(i-2)/2/h;endu0 = u1; end u = u1(2*M+1):(2*M+n); format short; 用 Richtmyer 多步格式解对流方程function u = peRich(a,dt,n,minx,maxx,M) format long; h = (maxx-minx)/(n-1);for j=1:(n+4*M)u0(j) = IniU(minx+(j-2*M-1)*h); end u1 = u0; for k=1:Mfor i=2*k+1:n+4*M-2*ktmpU1 = -dt*a*(u0(i+2)-u0(i)/h/4+(u0(i+2)+u0(i)/2;tmpU2 = -dt*a*(u0(i)-u0(i-2)/h/4+(u0(i)+u0(i-2)/2; u1(i) = -dt*a*(tmpU1-tmpU2)/h/2+u0(i);endu0 = u1; end u = u1(2*M+1):(2*M+n); format short; 用拉克斯-温德洛夫多步格式解对流方程function u = peMLW(a,dt,n,minx,maxx,M) format long; h = (maxx-minx)/(n-1); for j=1:(n+2*M)u0(j) = IniU(minx+(j-M-1)*h); end u1 = u0; for k=1:Mfor i=k+1:n+2*M-ktmpU1 = -dt*a*(u0(i+1)-u0(i)/h/2+(u0(i+1)+u0(i)/2;tmpU2 = -dt*a*(u0(i)-u0(i-1)/h/2+(u0(i)+u0(i-1)/2; u1(i) = -dt*a*(tmpU1-tmpU2)/h+u0(i);endu0 = u1; end u = u1(M+1):(M+n); format short; 用 MacCormack 多步格式解对流方程function u = peMC(a,dt,n,minx,maxx,M) format long; h = (maxx-minx)/(n-1); for j=1:(n+2*M)u0(j) = IniU(minx+(j-M-1)*h); end u1 = u0; for k=1:Mfor i=k+1:n+2*M-ktmpU1 = -dt*a*(u0(i+1)-u0(i)/h+u0(i);tmpU2 = -dt*a*(u0(i)-u0(i-1)/h+u0(i-1); u1(i) = -dt*a*(tmpU1-tmpU2)/h/2+(u0(i)+tmpU1)/2;endu0 = u1; end u = u1(M+1):(M+n); format short; 用拉克斯用拉克斯-弗里德里希斯格式解二维对流方程的初值问题弗里德里希斯格式解二维对流方程的初值问题function u = pe2LF(a,b,dt,nx,minx,maxx,ny,miny,maxy,M) %啦-佛format long; hx = (maxx-minx)/(nx-1); hy = (maxy-miny)/(ny-1); for i=1:nx+2*Mfor j=1:(ny+2*M)u0(i,j) = Ini2U(minx+(i-M-1)*hx,miny+(j-M-1)*hy);end end u1 = u0; for k=1:Mfor i=k+1:nx+2*M-kfor j=k+1:ny+2*M-ku1(i,j) = (u0(i+1,j)+u0(i-1,j)+u0(i,j+1)+u0(i,j-1)/4 .-a*dt*(u0(i+1,j)-u0(i-1,j)/2/hx .-b*dt*(u0(i,j+1)-u0(i,j-1)/2/hy;endendu0 = u1; end u = u1(M+1):(M+nx),(M+1):(M+ny); format short; 用拉克斯用拉克斯-弗里德里希斯格式解二维对流方程的初值问题弗里德里希斯格式解二维对流方程的初值问题function u = pe2FL(a,b,dt,nx,minx,maxx,ny,miny,maxy,M) %近似分裂format long; hx = (maxx-minx)/(nx-1); hy = (maxy-miny)/(ny-1); for i=1:nx+4*Mfor j=1:(ny+4*M)u0(i,j) = Ini2U(minx+(i-2*M-1)*hx,miny+(j-2*M-1)*hy);end end u1 = u0;for k=1:Mfor i=2*k+1:nx+4*M-2*kfor j=2*k-1:ny+4*M-2*k+2tmpU(i,j) = u0(i,j) - a*dt*(u0(i+1,j)-u0(i-1,j)/2/hx + .(a*dt/hx)2*(u0(i+1,j)-2*u0(i,j)+u0(i-1,j)/2;endendfor i=2*k+1:nx+4*M-2*kfor j=2*k+1:nx+4*M-2*ku1(i,j) = tmpU(i,j) - b*dt*(tmpU(i,j+1)-tmpU(i,j-1)/2/hy + .(b*dt/hy)2*(tmpU(i,j+1)-2*tmpU(i,j)+tmpU(i,j-1)/2;endendu0 = u1; end u = u1(2*M+1):(2*M+nx),(2*M+1):(2*M+ny); format short;
收藏 下载该资源
网站客服QQ:2055934822
金锄头文库版权所有
经营许可证:蜀ICP备13022795号 | 川公网安备 51140202000112号