FEM - lax friedrichs two step method

[blizzard12345]

hi i have written a code which according to my tutor is the lax friedrichs two step technique , however i can't see how to show how i came across this code (i kinda changed random things and hoped at some point he would say its correct.

thanks for any help in advance kyle :)

source for lax Friedrich method: http://www.scribd.com/doc/49845422/33/Lax-Wendroff-Method (page 69, i understand that it is only showing the one step for this source)

my code: (the program in question is a MATLAB script)



while time <= output
%calculate stability CFL condition to ensure stability
%(dt/dx)*(velocity + root(2*g))
wmax=0;
for i = 1:nx
w=abs(uh(i)/h(i)) + sqrt(h(i)*abs(g));
if(wmax<w)wmax=w;
end
end
dt=cfl*dx/wmax;
time = time + dt;
time
%calculate half
for i = 1:nx-1
h_half(i+1) = 0.5*((couple_one(h(i+1),uh(i+1)) + couple_one(h(i),uh(i))) - ((dx/dt)*( h(i+1) - h(i) )));
uh_half(i+1) = 0.5*((couple_two(h(i+1),uh(i+1),g) + couple_two(h(i),uh(i),g)) - ((dx/dt)*( uh(i+1) - uh(i) )));


end
%full time step calcs
for i = 2:nx-1


h(i) = h(i) - (dt/dx)*(h_half(i+1) - h_half(i));
uh(i) = uh(i) - (dt/dx)*(uh_half(i+1) - uh_half(i));
end
%update big H ready for plot
for i = 1:nx
H(i) = h(i);
UH(i) = uh(i);
end

 

[gtoguy97]

plot(x,H,'r'); hold on endThis code is based on the Lax Friedrichs two step technique. The first step is a calculation of the stability CFL condition to ensure that the numerical method is stable. This is done by calculating the maximum wave speed in the domain and setting the time step size, dt, according to the CFL condition (dt/dx)*(velocity + root(2*g)). The second step is to calculate the half step values for the height and velocity in the domain. This is done by using the couple_one and couple_two functions to calculate the values at each point. The third step is to calculate the full time step values for the height and velocity in the domain using the half step values from the previous step. This is done by subtracting (dt/dx)*(h_half(i+1) - h_half(i)) from the height and (dt/dx)*(uh_half(i+1) - uh_half(i)) from the velocity. Finally, the fourth step is to update the height and velocity values in the domain and plot them. The values are updated by setting the H and UH vectors equal to the h and uh vectors and then plotting the values in the H vector. Overall, this code follows the steps of the Lax Friedrichs two step technique for solving hyperbolic partial differential equations numerically.