对流方程各种格式代码matlab

对流方程各种格式代码matlab

对流方程——偏微分方程的数值解法

用迎风格式解对流方程

function u = peYF(a,dt,n,minx,maxx,M)

format long;

h = (maxx-minx)/(n-1);

if a>0

for 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:M

if a>0

for i=(k+1):n+M

u1(i) = -dt*a*(u0(i)-u0(i-1))/h+u0(i);

end

else

for i=1:n+M-k

u1(i) = -dt*a*(u0(i+1)-u0(i))/h+u0(i);

end

end

u0 = u1;

end

if a>0

u = u1((M+1):M+n);

else

u = 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:M

for i=k+1:n+2*M-k

对流方程各种格式代码matlab

u1(i) = -dt*a*(u0(i+1)-u0(i-1))/h/2+(u0(i+1)+u0(i-1))/2; end

u0 = 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:M

for i=k+1:n+2*M-k

u1(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);

end

u0 = 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:M

for i=2*k+1:n+2*M

u1(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;

end

u0 = 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);

对流方程各种格式代码matlab

for j=1:(n+4*M)

u0(j) = IniU(minx+(j-2*M-1)*h);

end

u1 = u0;

for k=1:M

for i=2*k+1:n+4*M-2*k

tmpU1 = -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);

end

u0 = 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:M

for i=k+1:n+2*M-k

tmpU1 = -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);

end

u0 = 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:M

for i=k+1:n+2*M-k

tmpU1 = -dt*a*(u0(i+1)-u0(i))/h+u0(i);

tmpU2 = -dt*a*(u0(i)-u0(i-1))/h+u0(i-1);

对流方程各种格式代码matlab

u1(i) = -dt*a*(tmpU1-tmpU2)/h/2+(u0(i)+tmpU1)/2; end

u0 = 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*M

for 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:M

for i=k+1:nx+2*M-k

for j=k+1:ny+2*M-k

u1(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;

end

end

u0 = 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*M

for 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;

对流方程各种格式代码matlab

for k=1:M

for i=2*k+1:nx+4*M-2*k

for j=2*k-1:ny+4*M-2*k+2

tmpU(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;

end

end

for i=2*k+1:nx+4*M-2*k

for j=2*k+1:nx+4*M-2*k

u1(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; end

end

u0 = u1;

end

u = u1((2*M+1):(2*M+nx),(2*M+1):(2*M+ny));

format short;