- #1
blizzard12345
- 14
- 0
i am trying to spproximate a PDE in the form below using the lax wendroff 2 step method in MATLAB coding:
[h ; hu ] = [ hu ; hu^2 + 1/2gh^2] = [0; -ghbz] (where bz will equal zero)
i believe this is then the case
d(h)/dt + d(hu)/dx = 0
and
d(hu)/dt + d(hu^2 + 1/2gh^2)/dx = 0
as the st venant equation was formed by depth deriving the 1D continuity equation
however i have attempted to program this and i keep produced a column matrix full of "NaN" could some one try to point out where i have been going wrong? thanks kyle
code for matlab:
%%%%%%%%%%%%%%%% function compare close all clc clear all
%intial values ntime = 1000000; dt=0.00050; nx=100; time = 0; a=2;output=0.24; %step size calculation dx= (1/nx); %create size of u_int vector u_int = zeros(nx,2); %create little u_int vector for the initial values of height beginning and %end depending which value of A is used u_int(nx,2) =1; u_int(1,1) = 1;
%loop for the two directions needed for vec = 1:2 % %lax_wendroff [H] = lax_wendroff_2(nx,a,output, dt, ntime, time,u_int, vec,dx); figure(vec) H plot(H(:,vec),'r*');grid on;legend('centered'); hold on end end
function [H] = lax_wendroff_2(nx,a,output, dt, ntime, time,u_int, vec,dx); %compute original size matracies g = 9.81; uh=zeros(nx,2); update = uh ; H = uh; Hy = uh; UH = uh; if vec == 1; %for i = 1:nx u_int(nx-2:nx,2) = 1; u_int(1:3,1) = 1; %end for i = 1:nx UH(i,vec)=u_int(i,vec); Hy(i,vec)=u_int(i,vec); uh(i,vec)=u_int(i,vec); H(i,vec)=u_int(i,vec); A=a; end else uh(nx,vec)=u_int(nx,vec); h(nx,vec)=u_int(nx,vec); A=-a; end for p = 1:ntime; if(time + dt>output); dt=output-time; end % calculate coefficent coeff=dt/dx; %calculate interim steps for i = 1:(nx-1) % first derivative Hy(i+1,vec)= ((H(i+1,vec)+H(i,vec))/2) - ((coeff*A/2)*((uh(i+1,vec) - uh(i,vec))));
UH(i+1,vec)=((uh(i,vec)+uh(i+1,vec))/2) ;
-coeff*A*(((uh(i+1,vec)^2)/(H(i+1,vec)+g/2*H(i+1,vec)^2)) - ((uh(i,vec)^2)/(H(i,vec)+g/2*H(i,vec)^2)));
end
%calculate full time steps
for i = 2:(nx-1)
update_h(i+1,vec) = H(i,vec) - coeff*A*(UH(i+1,vec) - UH(i,vec));%%%% check i notation and convention may be wrong way
update_uh(i+1,vec) = (uh(i,vec)) - (coeff*A*(((UH(i+1,vec)^2)/(Hy(i+1,vec)+g/2*Hy(i+1,vec)^2))- ((UH(i,vec)^2)/(Hy(i,vec)+g/2*Hy(i,vec)^2))));
end
for i = 3:(nx-1)
H(i,vec) = update_h(i,vec);
uh(i,vec) = update_uh(i,vec);
end
time = time + dt;
if time==output
break
end
end end %%%%%%%%%%%%%%%%%
[h ; hu ] = [ hu ; hu^2 + 1/2gh^2] = [0; -ghbz] (where bz will equal zero)
i believe this is then the case
d(h)/dt + d(hu)/dx = 0
and
d(hu)/dt + d(hu^2 + 1/2gh^2)/dx = 0
as the st venant equation was formed by depth deriving the 1D continuity equation
however i have attempted to program this and i keep produced a column matrix full of "NaN" could some one try to point out where i have been going wrong? thanks kyle
code for matlab:
%%%%%%%%%%%%%%%% function compare close all clc clear all
%intial values ntime = 1000000; dt=0.00050; nx=100; time = 0; a=2;output=0.24; %step size calculation dx= (1/nx); %create size of u_int vector u_int = zeros(nx,2); %create little u_int vector for the initial values of height beginning and %end depending which value of A is used u_int(nx,2) =1; u_int(1,1) = 1;
%loop for the two directions needed for vec = 1:2 % %lax_wendroff [H] = lax_wendroff_2(nx,a,output, dt, ntime, time,u_int, vec,dx); figure(vec) H plot(H(:,vec),'r*');grid on;legend('centered'); hold on end end
function [H] = lax_wendroff_2(nx,a,output, dt, ntime, time,u_int, vec,dx); %compute original size matracies g = 9.81; uh=zeros(nx,2); update = uh ; H = uh; Hy = uh; UH = uh; if vec == 1; %for i = 1:nx u_int(nx-2:nx,2) = 1; u_int(1:3,1) = 1; %end for i = 1:nx UH(i,vec)=u_int(i,vec); Hy(i,vec)=u_int(i,vec); uh(i,vec)=u_int(i,vec); H(i,vec)=u_int(i,vec); A=a; end else uh(nx,vec)=u_int(nx,vec); h(nx,vec)=u_int(nx,vec); A=-a; end for p = 1:ntime; if(time + dt>output); dt=output-time; end % calculate coefficent coeff=dt/dx; %calculate interim steps for i = 1:(nx-1) % first derivative Hy(i+1,vec)= ((H(i+1,vec)+H(i,vec))/2) - ((coeff*A/2)*((uh(i+1,vec) - uh(i,vec))));
UH(i+1,vec)=((uh(i,vec)+uh(i+1,vec))/2) ;
-coeff*A*(((uh(i+1,vec)^2)/(H(i+1,vec)+g/2*H(i+1,vec)^2)) - ((uh(i,vec)^2)/(H(i,vec)+g/2*H(i,vec)^2)));
end
%calculate full time steps
for i = 2:(nx-1)
update_h(i+1,vec) = H(i,vec) - coeff*A*(UH(i+1,vec) - UH(i,vec));%%%% check i notation and convention may be wrong way
update_uh(i+1,vec) = (uh(i,vec)) - (coeff*A*(((UH(i+1,vec)^2)/(Hy(i+1,vec)+g/2*Hy(i+1,vec)^2))- ((UH(i,vec)^2)/(Hy(i,vec)+g/2*Hy(i,vec)^2))));
end
for i = 3:(nx-1)
H(i,vec) = update_h(i,vec);
uh(i,vec) = update_uh(i,vec);
end
time = time + dt;
if time==output
break
end
end end %%%%%%%%%%%%%%%%%