偏微分差分四种格式的matlab程序.doc

偏微分差分四种格式的matlab程序

(1)Lax-Fridrichs格式

function u = LaxFridrichs(a,dt,n,x1,x2,m)

format long;

h = (x2-x1)/(n-1);

for j=1:(n+2*m)

u0(j) =U(x1+(j-m-1)*h);

end

u1 = u0;

for k=1:m

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

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

end

u0 = u1;

end

x=0.01:0.01:n/100;

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

plot(x,u,'r');

axis([0 1 -0.5 1.5]);

hold on

y=1.*(x<=0.5)+0.*(x>0.5);

plot(x,y,'c')

title('Lax-Fridrichs')

format short;

function ux =U(x)

format long;

if x<=0

ux=1;

else

ux=0;

end

>> u = LaxFridrichs(1,0.005,101,0,1,100)

运行结果

(2)迎风格式

function u = Yingfeng(a,dt,n,x1,x2,m)

format long;

h = (x2-x1)/(n-1);

if a>0

for j=1:(n+m)

u0(j) =U(x1+(j-m-1)*h);

end

else

for j=1:(n+m)

u0(j) =U(x1+(j-1)*h);

end

end

u1 = u0;

for k=1:m

if a>0

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

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

end

else

for i=1:n+m-k

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

end

end

u0 = u1;

end

x=0.01:0.01:n/100;

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

plot(x,u,'r');

axis([0 1 -0.5 1.5]);

hold on

y=1.*(x<=0.5)+0.*(x>0.5);

plot(x,y,'c')

format short;

title('迎风格式')

function ux =U(x)

format long;

if x<=0

ux=1;

else

ux=0;

end

>> u=Yingfeng(1,0.005,101,0,1,100)

运行结果

(3)Lax-Wendroff格式

function u = LaxWendroff(a,dt,n,x1,x2,m)

format long;

h = (x2-x1)/(n-1);

for j=1:(n+2*m)

u0(j) =U(x1+(j-m-1)*h);

end

u1 = u0;

for k=1:m

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

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

end

u0 = u1;

end

x=0.01:0.01:n/100;

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

plot(x,u,'r');

axis([0 1 -0.5 1.5]);

hold on

y=1.*(x<=0.5)+0.*(x>0.5);

plot(x,y,'c')

title('Lax-Wendroff')

format short

function ux =U(x)

format long;

if x<=0

ux=1;

else

ux=0;

end

>> u= LaxWendroff(1,0.005,101,0,1,100)

运行结果

(4)Beam-Warming格式

function u =BeamWarming(a,dt,n,x1,x2,m)

format long;

h = (x2-x1)/(n-1);

for j=1:(n+2*m)