Error in graph - Finite difference code

I'm a research scholar and I've started coding in MATLAB recently. I'm trying to solve a system of differential equations using finite difference method. The paper I'm working on is by Sharma et al. (I've given the link to the paper here Sharma et al. (2019)). In the paper, they have solved the finite difference equations using Thomas algorithm. But, I'm trying to solve the same finite difference equations without using Thomas algorithm. I've attached my code for your reference. There are no errors shown in the code. However, I'm unable to recreate the graphs given in the paper. May I know where I've gone wrong?

Nr=21; Nz=161;                                  %number of values for radius r, z
R=1; L=16;                                      %max values along each axis
dr=R/(Nr-1); dz=L/(Nz-1);                       %deltas
r=(0:Nr-1)*dr; z=(0:Nz-1)*dz;                   %vectors for r, z
u=zeros(Nr,Nz);  w=zeros(Nr,Nz); 
p=zeros(Nr,Nz);  theta=zeros(Nr,Nz); 
sigma=zeros(Nr,Nz);                             %allocate arrays

R=1; M=1; delta=0.6; epsilon=0.005;             %constants
alpha=0; n=2; beta=1; Df=1; Sc=0.5; Sr=0.2;     %constants


eta=(delta*(n^(n/(n-1))))/(n-1);
%alpha=a/b
h=(1+(epsilon*beta))*(1-eta*((beta-alpha)-(beta-alpha)^n));

%w(:,1)=w*ones(Nr,1);     %w distribution at z=0
%w(1,:)=w*ones(1,Nz);     %w distribution at r=0


    for j=1:Nz-1
        for i=2:h
            u(i+1,j)=(1/dz*(r(i)))*(((u(i,j))*(r(i))*dz)-((u(i,j))*(dr)*(dz))+((w(i,j))*(r(i))*(dr))-((w(i,j+1))*(r(i))*(dr)));
            %(u(i+1,j)-u(i,j))/(dr)+(u(i,j)/r(i,j))+(w(i,j+1)-w(i,j))/(dz)==0;
            p(i+1,j)=p(i,j);
            p(i,j+1)=p(i,j)+(dz/r(i))*((w(i+1,j)-w(i,j))/dr)+((dz)/(dr)^2)*(w(i+1,j)-2*w(i,j)+w(i-1,j))-(M^2)*dz*w(i,j);
            ((1+(4/3*R))*(((theta(i+1,j)-theta(i,j))/(dr*(r(i))))+((theta(i+1,j)-2*theta(i,j)+(theta(i-1,j)))/((dr)^2))))+(Df*(((sigma(i+1,j)-sigma(i,j))/(dr*(r(i))))+((sigma(i+1,j)-2*sigma(i,j)+sigma(i-1,j))/((dr)^2))))==0;
            ((1/Sc)*(((sigma(i+1,j)-sigma(i,j))/(dr*(r(i))))+((sigma(i+1,j)-2*sigma(i,j)+sigma(i-1,j))/((dr)^2))))+((Sr)*(((theta(i+1,j)-theta(i,j))/(dr*(r(i))))+((theta(i+1,j)-2*theta(i,j)+theta(i-1,j))/((dr)^2))))==0;
        
        w(1,j)=w(2,j); theta(1,j)=theta(2,j); sigma(1,j)=sigma(2,j);  %boundary condition at r=0
        w(Nr,j)=0; theta(Nr,j)=0; sigma(Nr,j)=0;   %boundary condition at r=R
        %w(i,1)=0; theta(i,1)=0; sigma(i,1)=0; 
        %w(i,Nz)=0; theta(i,Nz)=0; sigma(i,Nz)=0;
        end
    end

%plot(r,[w(10,:);w(15,:);w(20,:)])
plot(r,theta)

[Subha S]

Hello all. I'm still trying to find what mistake I've made in the code above. Is the error due to no boundary conditions along z direction (along j direction)? Is that why I'm getting the graph as a line at the axis and not any other curve? Or have I made any other mistake in the code? Can anyone please help me find where I went wrong?

[Torsten]

Do you understand which variables are solved for in the article ? I see w, theta, p and sigma as discretized variables, but only two equations on page 330.

[Subha S]

Hello Torsten!!

I now found from another paper that dp/dz is taken as constant =-P0=-10. This means that I'll have 3 equations in three unknowns - w, theta and sigma. Further, since there are no more z terms/ terms that depend on z, there is no need for j loop. Only i loop is needed. However, my problem persits as earlier. I get the same graph as before after making the corrections to my code. I'm attaching my code below. I'm not understanding where I'm going wrong.

Nr=1001; Nz=1001;                               %number of values for radius r, z
S=1; L=1;                                       %max values along each axis
dr=S/(Nr-1); dz=L/(Nz-1);                       %deltas
r=(0:Nr-1)*dr; z=(0:Nz-1)*dz;                   %vectors for r, z
%u=zeros(Nr,Nz);  
w=zeros(Nr); 
%p=zeros(Nr,Nz);  
theta=zeros(Nr); 
sigma=zeros(Nr);                             %allocate arrays

R=1; M=1; delta=0.6; epsilon=0.005; P_0 =10;    %constants
alpha=0; n=2; beta=1; Df=1; Sc=0.5; Sr=0.2;     %constants


eta=(delta*(n^(n/(n-1))))/(n-1);
%alpha=a/b
h=(1+(epsilon*beta))*(1-eta*((beta-alpha)-(beta-alpha)^n));

%w(:)=w*ones(Nr);     %w distribution at z=0
%theta(:)=theta*ones(Nr); sigma(:)=sigma*ones(Nr);

    %for j=1:Nz-1
        for i=2:h
            %u(i+1,j)=(1/dz*(r(i)))*(((u(i,j))*(r(i))*dz)-((u(i,j))*(dr)*(dz))+((w(i,j))*(r(i))*(dr))-((w(i,j+1))*(r(i))*(dr)));
            %(u(i+1,j)-u(i,j))/(dr)+(u(i,j)/r(i,j))+(w(i,j+1)-w(i,j))/(dz)==0;
            %p(i+1,j)=p(i,j);
            -P_0==(1/r(i))*((w(i+1)-w(i))/dr)+((1)/(dr)^2)*(w(i+1)-2*w(i)+w(i-1))-(M^2)*w(i);
            ((1+(4/3*R))*(((theta(i+1)-theta(i))/(dr*(r(i))))+((theta(i+1)-2*theta(i)+(theta(i-1)))/((dr)^2))))+(Df*(((sigma(i+1)-sigma(i))/(dr*(r(i))))+((sigma(i+1)-2*sigma(i)+sigma(i-1))/((dr)^2))))==0;
            ((1/Sc)*(((sigma(i+1)-sigma(i))/(dr*(r(i))))+((sigma(i+1)-2*sigma(i)+sigma(i-1))/((dr)^2))))+((Sr)*(((theta(i+1)-theta(i))/(dr*(r(i))))+((theta(i+1)-2*theta(i)+theta(i-1))/((dr)^2))))==0;
        
        w(1)=w(2); theta(1)=theta(2); sigma(1)=sigma(2);  %boundary condition at r=0
        w(Nr)=0; theta(Nr)=0; sigma(Nr)=0;   %boundary condition at r=R
        
        %w(i,1)=w(i,2); theta(i,1)=theta(i,2); sigma(i,1)=sigma(i,2); 
        %w(i,Nz)=0; theta(i,Nz)=0; sigma(i,Nz)=0;
        end
    %end

%plot(r,[w(10,:);w(15,:);w(20,:)])
plot(r,theta)

[Torsten]

This means that I'll have 3 equations in three unknowns - w, theta and sigma

Why three equations ? I only see two equations for three unknowns w, theta and sigma.

And p is not constant in z-direction as clearly written in the first equation.

[Subha S]

Why three equations ? I only see two equations for three unknowns w, theta and sigma.

Sir, there are three equations. I'm attaching the screenshot of the equations from the paper. The first equation is =..., second equation is and the third equation is .

As for pressure, I've found that most papers have considered pressure along the z direction to be a constant, say, . And the continuity equation for pressure is not considered for calculation since the pressure along r-direction is negligible. This is why finite difference is not applied for the term . This explanation is not given in this paper but I found the explanation from other papers.

Therefore, I need to solve for w, theta and sigma along r direction alone.

[Torsten]

Ok, now I see that the discretized algebraic equations refer to equations (10), (11) and (12) of the paper.

But I cannot believe that P along the z-direction is considered constant. If this were the case, there would be no flow neither in r nor in z direction. Maybe you mean that dP/dz is considered constant.

Although your equations are linear in the unknowns, it will make a lot of effort to collect terms such that your system reads A*[w ;sigma;theta] = b and that you could solve for the unknown vector [w; sigma; theta]. So I suggest that you just implement the system as written in the paper (f(w,sigma,theta) = 0) and use "fsolve" to solve for the unknowns. The (nonlinear) solver "fsolve" should converge in one step to the solution because your system is linear.

Or if you have problems with the discretization, simply use "bvp4c" or "bvp5c" for the solution of stationary boundary value problems.

[Subha S]

Thank you for your reply Sir.

I'll try with fsolve or bvp4c/bvp5c and get back to you.

[Torsten]

Equations (11) and (12) are in principle a homogenous system of equations in the unknowns 1/r * d/dr (r*dtheta/dr) and 1/r * d/dr(r*dsigma/dr). The solution is

1/r * d/dr (r*dtheta/dr) = 0
1/r * d/dr (r*dsigma/dr) = 0 

So with the boundary conditions given, it's hard for me to understand why not theta = sigma = 0 for all r follows.

[Subha S]

You are correct in the regard that theta and sigma equations are essentially homogeneous system. But if theta and sigma are zeros for all r, then how did they obtain non-zero graphs for those variables?

[Torsten]

But if theta and sigma are zeros for all r, then how did they obtain non-zero graphs for those variables?

Good question - I don't know.

[Subha S]

Hi Sir. I'm now working on another paper. I'll see how I get the graphs for that paper. If I get struck even in this paper, I'll ask question in this forum. Thank you Torsten Sir for spending time and commenting. Your comments are useful. Thank you once again

[Subha S]

Hello, Torsten.

Now, I'm solving the same set of equations using bvp4c solver as suggested by you. I've split the code into two codes - one for w equation and the other for theta and sigma. I'm getting the graphs, however, I'm not able to produce the graphs given in the paper. I tried to set axis limits, but I got error. How to get the desired graphs?

R=4;
P_0=10;
k=1;
eps=1.0000e-06;
r=linspace(0,R,10);
solinit = bvpinit(linspace(0,R,10),[1,1]);  
options = bvpset('RelTol', 10e-4,'AbsTol',10e-7);
sol = bvp4c(@(r,y)soret(r,y,P_0),@soret_BC,solinit,options);
y=deval(sol,r);
plot(r,y(1,:))

function dydx = soret(r,y,P_0) % Details ODE to be solved
        dydx = zeros(2,1);
        dydx(1) = y(2);
        dydx(2) = -(P_0)-((1/(r+eps))*y(2))+(y(1)) ; % This equation is invalid at r = 0, so taken as r+eps
end
    
function res = soret_BC(ya,yb) % Details boundary conditions
        res=[ya(2);yb(1)];
end

R=6;
Q=1;
D_f=1;
Sc=0.5;
Sr=0.2;
k=1;
eps=1.0000e-06;
r=linspace(0,R,100);
solinit = bvpinit(linspace(0,R,100),[1,1,1,1]); 
options = bvpset('RelTol', 10e-4,'AbsTol',10e-7);
sol = bvp4c(@(r,y)soret(r,y,Q,D_f,Sc,Sr,k),@soret_BC,solinit,options);
y=deval(sol,r);
plot(r,y(2,:))

function dydx = soret(r,y,Q,D_f,Sc,Sr,k) % Details ODE to be solved
        dydx = zeros(4,1);
        dydx(1) = y(3); dydx(2)=y(4);
        dydx(3) = (1)/(1+((4)/(3*Q)))*(((-(1/k))*(y(3))*(1+((4)/(3*Q))))+((-D_f)*((1/(r+eps)))*(y(4)))+((-D_f)*(dydx(4))));
        dydx(4) = Sc*(((-1/Sc)*(1/k)*y(4))+((-Sr)*(1/(r+eps))*(y(3)))+((-Sr)*(dydx(3))));
          % This equation is invalid at r = 0, so taking 1/r as 1/k and
          % 1/r+eps
end
    
function res = soret_BC(ya,yb) % Details boundary conditions
        res=[ya(3); ya(4); yb(1); yb(2)];
end

Function definitions in a script must appear at the end of the file.
Move all statements after the "soret_BC" function definition to before the first local function definition.

[Torsten]

Which equations do you use ?

For the first equation, I get

d^2w/dr^2 + 1/r * dw/dr - M^2*w = dPdz,

thus

dy(1) = y(2);
dy(2) = -1/r * y(2) + M^2 * y(1) + dPdz

I cannot recover that from your code.

Note that you don't need to care about r = 0 in the equations setting if you exclude r = 0 in the mesh definition:

r=linspace(eps,R,10);

[Subha S]

I've attached the equations that I have used in the bvp4c code.

Note that you don't need to care about r = 0 in the equations setting if you exclude r = 0 in the mesh definition

Thank you for pointing this !!

I hope that the equations I used are correct. But I don't know how to get the same graphs as given in the paper.

[Torsten]

As I already wrote, the equations (2) and (3) from the paper (and thus yours in the bvp code) must be wrong.

I showed you that they reduce to

1/r * d/dr (r*dtheta/dr) = 0
1/r * d/dr (r*dsigma/dr) = 0 

and that they give theta = sigma = 0 (together with the boundary conditions given).

bvp4c doesn't work for combinations of second-order derivatives. You explicitly have to solve for the second-order derivatives and then reduce the system to a first-order system.

For your equations, this means that they can only be solved in the form

1/r * d/dr (r*dtheta/dr) = 0
1/r * d/dr (r*dsigma/dr) = 0 

not as

a * 1/r * d/dr (r*dtheta/dr) + b * 1/r * d/dr (r*dsigma/dr) = 0
c * 1/r * d/dr (r*dtheta/dr) + d * 1/r * d/dr (r*dsigma/dr) = 0

[Subha S]

Thank you Sir for pointing out the error. I've changed the equations for theta and sigma. However, I'm not obtaining a smooth curve, there is a discontinuity that arises due to the problem itself. How to rectify the discontinuity? The same thing happens with w graph as well. Please suggest on how to proceed. I'm attaching my code for both w and theta, sigma.

R=1.5;
P_0=10;
k=1;
eps=1.0000e-06;
r=linspace(-R,R,1000);
solinit = bvpinit(linspace(-R,R,10),[0,0]);  
options = bvpset('RelTol', 10e-4,'AbsTol',10e-7);
sol = bvp4c(@(r,y)soret(r,y,P_0),@soret_BC,solinit,options);
y=deval(sol,r);
plot(r,y(1,:))

function dydx = soret(r,y,P_0) % Details ODE to be solved
        dydx = zeros(2,1);
        dydx(1) = y(2);
        dydx(2) = -(P_0)-((1/(r+eps))*y(2))+(y(1)) ; % This equation is invalid at r = 0, so taken as r+eps
end
    
function res = soret_BC(ya,yb) % Details boundary conditions
        res=[ya(2);yb(1)];
end

The code for theta and sigma are given below

R=1.5;

Function definitions in a script must appear at the end of the file.
Move all statements after the "soret_BC" function definition to before the first local function definition.

Q=1;
D_f=1;
Sc=0.5;
Sr=0.2;
k=1;
eps=1.0000e-06;
r=linspace(-R,R,100);
solinit = bvpinit(linspace(-R,R,100),[0,0,1,1]); 
options = bvpset('RelTol', 10e-4,'AbsTol',10e-7);
sol = bvp4c(@(r,y)soret(r,y),@soret_BC,solinit,options);
y=deval(sol,r);
plot(r,y(2,:))

function dydx = soret(r,y) % Details ODE to be solved
        dydx = zeros(4,1);
        dydx(1) = y(3); dydx(2)=y(4);
        dydx(3) = -(y(3))/((r+eps));
        dydx(4) = -(y(4))/((r+eps));
          % This equation is invalid at r = 0
end
    
function res = soret_BC(ya,yb) % Details boundary conditions
        res=[ya(3); ya(4); yb(1); yb(2)];
end

[Torsten]

R=1.5;
P_0=10;
k=1;
eps=1.0000e-04;
r=linspace(eps,R,1000);
solinit = bvpinit(r,[0,0]);  
options = bvpset('RelTol', 10e-4,'AbsTol',10e-7);
sol = bvp4c(@(r,y)soret(r,y,P_0),@soret_BC,solinit,options);
y=deval(sol,r);
plot(r,y(1,:))

function dydx = soret(r,y,P_0) % Details ODE to be solved
        dydx = zeros(2,1);
        dydx(1) = y(2);
        dydx(2) = -P_0 - 1/r * y(2) + y(1) ; % This equation is invalid at r = 0, so taken as r+eps
end
    
function res = soret_BC(ya,yb) % Details boundary conditions
        res=[ya(2);yb(1)];
end



R=1.5;
Q=1;
D_f=1;
Sc=0.5;
Sr=0.2;
k=1;
eps=1.0000e-04;
r=linspace(eps,R,100);
solinit = bvpinit(r,[0,0,1,1]); 
options = bvpset('RelTol', 1e-8,'AbsTol',1e-8);
sol = bvp4c(@(r,y)soret(r,y),@soret_BC,solinit,options);
y=deval(sol,r);
plot(r,y(2,:))

function dydx = soret(r,y) % Details ODE to be solved
        dydx = zeros(4,1);
        dydx(1) = y(3); dydx(2)=y(4);
        dydx(3) = -y(3)/r;
        dydx(4) = -y(4)/r;
          % This equation is invalid at r = 0
end
    
function res = soret_BC(ya,yb) % Details boundary conditions
        res=[ya(3); ya(4); yb(1); yb(2)];
end