1-D convection equation
As we are trying to solve the 1-D wave/convection equation, we need to know how we are solving and what kind of procedures we're foloowing to get the solution.
First of all the 1-D convective equation falls under the category of hyperbolioc partial differential equation. To solve hyperbolic equation we can follow different kinds of numerical schemes. In this particular hyperbolic equation, we can solve this by three different methods. FTCS, FTBS and FTFS. But here we're following FTBS, it stands for Forward Time Backward Space. In other words, we discretize the equation usinf forward scheme for time derivative and backward scheme for space derivative term.
clear all
close all
clc
%inputs for linear convection
L=1;
c=1;
n=80; % number of nodes
%initialize arrays
x=linspace(0,L,n);
u=ones(1,n);
## dx=L/(n-1)
dx=x(2)-x(1);
%calculations and defining start and end posiions for step function
for (j=1:20)
if (x(j)>=0.1 && x(j)<=0.3)
u(j)=2;
else
u(j)=1;
endif
endfor
%BC
uold=u;
u_initial=u;
grid on
nt=40;
t_t=0.4;
ct=1;
for k=1:nt %time loop
dt=(t_t/nt);
% space loop
for i=2:n
u(i)=uold(i)-[(c*dt/dx)*(uold(i)-uold(i-1))];
endfor
% update old values
uold=u;
plot(x,u,'g')
axis([0 1 0 2])
xlabel('x_ length of domain')
ylabel('velocity')
pause(0.001)
text=sprintf('linear_convection_%d.png',ct);
title('Linear convection with 160 nodes');
saveas(1, text);
ct=ct+1;
endfor
% plot for comparing initial and final values
figure 2
plot(x, u_initial,'r')
hold on
plot(x, u,'-@')
axis([0 1 0 2])
xlabel('x_ length of domain');
ylabel('velocity');
legend('initial','final')
In the code provided above we can see that, for solving a partial differential euation of this kind we need two kinds of conditons, first one is boundary conditions and the other is initial condition. As we start with the code, we are assinging the values of length of domain, number of nodes this length is going to be divided and the total time of which it is going to be performed. After giving those values, we're giving function and values of 'u' which is velocity respective to the length of domain.
As we see the number plot in the above figures it showing the convection for 1 metre of length of 0.4 seconds long with 20 nodes so the space step is going to be 1/19=0.05263 and time step is 0.01 with 40 cycles. But we can observe that the velocity is not reached upto 2 m/s at any point of time.
If we compare the plot that shows initial and final velocities, we can observe that as we increase the number of nodes the funtion is becoming more sharp and steeper, meaning that we're going in the right direction of the condition what we provided in the code.
Now we can clearly say that as we increase the number of nodes then the function what we want is becoming more sharp and steeper. In this th velocity can even touch the peak value till this at 1.97 m/s.
As we saw in the previous plots, if we increase number of nodes the plot becoming more steeper and sharp. But contrast to the plots whatever we saw till now, here we can clearly see that the solution is blowing up and we're not getting what we desired. This behaviour is called instability of a scheme. So why is this happening? What is the reason behind that this solution is blowing up ?
For a hyperbolic equation to be stable and convergent we have a particular criteria called CFL number. As it stands for Courant-Fredrich-Lewis criteria.
`CFL=(c*(deltat)/(deltax))` here the c stands for courant number. It has a significance in determining the stability of the equation. Courant number physically tells you the speed at which the disturbance is travelling in a wave. When it is multiplied with the time step it will give you the distance of the disturbance in the travelling wave. So the CFL criteria is only valid for linear partial differntial equations. For numerical scheme to be stable this CFL number has to be less than or equal to 1. If the number is less than 1 then it shows that it is convergent as well as stable. If it is equal to 1 then it says that the equation is stable without any convergance or error.
As we discussed above the courant nnumber says at which the disturbance is propagating thtough the wave, and if it is multiplied with time step it shows how much distance the disturbance travelled in the wave. If we give a grid size that is too small for a wave to propogate in the time, then the results will blow up and we'll get the results whicha are unwanted with huge maginification. This happened in the fourth plot that shows the entire conevtion with 160 number of nodes, meaning that the CFL criteria is greater than 1 for 160 number of nodes. If we calculate the grid size it is nearly 0.006289, for this smaller grid size the wave cannot propagate, that is why we're getting the unwanted results ans blowned up plots. In this we can observe that as we increase the grid size, we always thought we are increasing the accuracy of the numerical scheme that are using, contrast to that we saw the blown up results ansd unwanted plots that are unstable. So, now we can conveniently say that grid size is always not the factor that decides the accuracy of the numerical scheme. Anything more can be too much, as in the case of grid size also.
So, we observed blown up results and convergent results, but where is the result that is neither decaying not increasing? Here we can give the number nodes as 100 and we can get the most valid result that we're looking for, this is because if we provide 1000 nodes then the CFL is exaclty equals 1 and hence the result.
As we know there are two other numerical schemes that we can apply to solve this equation, those are FTFS and FTCS. These two schemes are unconditonally unstable because the amplification factor is greater than 1 always. In simple words amplification factor is the ratio of error at n+1 th position or grid to the n th position or grid. This error is generated becuse there is difference betweeen approximate difference equation and exact difference equation. These two schemes FTFS and FTCS will always yield an amplification factor greater than 1 without any condtion.
Here are the simulations for both FTFS and FTCS schems