Grid Covergence

Grid Covergence

Solution 

clc

clear all

rho=1.2

At=0.1

L=1

nx=10

nx1=nx+1

hx=L/nx

xvec=-(1-hx*(0:nx))

alpha = At*sqrt(1+(xvec./L).^2)

beta=(At*xvec/(L*L))./sqrt(1+(xvec./L).^2)

p= beta./alpha

p1=p.*hx

D(1,1)=-1

D(1,2)=1

D(nx1-1,nx1-2)=1

D(nx1-1,nx1-1)=-(2+p1(nx1-1))

for i=2:nx-1

D(i,i-1)=1

D(i,i)=-(2+p1(i))

D(i,i+1)=(1+p1(i))

end

B=zeros(nx,1)

B(1,1)=hx

%B=[hx 0 0 0]’

X=inv(D)*B

X(nx1)=0 % final Bc as phi

% X denotes phi

u(1,1)=1;   % initial BC                                                                                                                                  for k=1:nx

u(k+1)=(X(k+1)-X(k))/hx

end

plot(xvec,u,’-b’,’LineWidth’,2,’DisplayName’,’ U ‘)

hold on

plot(xvec,X’,’-r’,’LineWidth’,2,’DisplayName’,’ phi ‘)

legend(‘show’)

xlabel(‘X’)

ylabel(‘u,phi’)

hold on