Scaling parameters in central difference solution

[robby991]

Hi, I developed Matlab code to solve the diffusion equation using the central difference equation, with an added term at the end. The equation is the following:

 dS/dt=Ds*d^2S/dx^2-(Vmax*S/Km+S)

In my code, the length of the space domain is very small, 1E-6. I would like to scale my code to run in the space domain between x = 0 and 1, rather than x = 0 to 1E-6. I propose to introduce a new variable, x2 = x/L. My question is what parameters/variables in my code need to be scaled also. My code is:

clear all;

numx = 10;                      %number of grid points in space
numt = 1000;                    %number of time steps to be iterated over 
tmax = .0045;
Length = 1E-6;                  %length of grid
Ds = .019E-9;                   %requirement Ds(dt)/dx^2 < .5, cm^2/sec
Vmax = 275E-6;                  %mol/cm^2sec
Km = 3E-3;                      %mol/cm^3

x = linspace(0,Length,numx);    %vector of x values, to be used for plotting
t = linspace(0,tmax,numt)';     %vector of t values, to be used for plotting
S = zeros(numt,numx);           %initialize everything to zero

dx = x(2)-x(1);                 %Define grid spacing in time
dt = t(2)-t(1);                 %Define grid spacing in time

%specify initial conditions%

t(1) = 0;      %1st t position = 0

S(1,:) = cos(x*pi/(2*Length));   
S(:,numx) = 0;

S_exact = cos(x*pi/(2*Length))*exp(-(pi*sqrt(Ds)/(2*Length))^2*tmax);

%iterate central difference equation% 

for j=1:numt-1  
    
      
    %2nd Derivative Central Difference Iteration% 
     
    for i=2:numx-1
      S(j+1,i) = S(j,i) + (dt/dx^2)*Ds*(S(j,i+1) - 2*S(j,i) + S(j,i-1))-((Vmax*dt*S(j,i))/(Km+S(j,i))); 
    end
    
    S(j+1,1)=S(j,1)+dt*Ds*2*(S(j,2)-S(j,1))./dx.^2-((Vmax*dt*S(j,i))/(Km+S(j,i)));      %Neumann Boundary Condition
    
end
   
plot(x,S(numt,:));  
hold on
plot(x,S_exact,'r*')

error = max(S(numt,:)-S_exact)

To scale, the domain would now be:

x2=x/Length;

and now, dx = x2(2)-x2(1);

How do I change me code to accommodate this? How do my other variables change/scale? i.e. Km, Ds, Vmax etc. Thank you.

 

[AlephZero]

This is exactly the same as working in different units, for example changing your lengths from km to mm (as an example of units which happen to be a factor of 10^6 different)

If you change the values of ALL the constants in your code which have physical dimensions, you don't need to change anything else.

 

[robby991]

I have one last question regarding this. Do I have to scale the matrix vectors in the time domain also? I scaled all the parameters related to the space domain (x), but I also have the time vectors (y). If so, do I divide by the same constant I used to scale the space domain?

 

[robby991]

I ran into a little difficulty and was wondering if my scaling of parameters was correct in my code. Say I have a grid of length "Length", and instead of going from 0 to Length I want to scale from 0 to 1. In addition I have the following constants in my calculations:

Length = 1E-4;                  %cm
Ds = .019E-9;                   % cm^2/sec
Vmax = 275E-6;                  %mol/cm^2sec
Km = 3E-3;                      %mol/cm^3

To scale from 0 to 1, I would simply make a new variable L where:

L = Length/Length.

Next, to scale the other variables I will do the following:

D = Ds/Length^2
V = Vmax*Length^2
K = Km*Length^3

Is this correct? I would appreciate any input on the matter. Thank you.

 

[blue_raver22]

I would recommend the following steps to scale your code:

1. Scale the space domain: As you have proposed, introduce a new variable x2 = x/L to scale the space domain from 0 to 1. This will require you to change the definition of dx as dx = x2(2)-x2(1).

2. Scale the diffusion coefficient (Ds): Since the diffusion coefficient is a function of space (Ds = Ds(x)), it will also need to be scaled accordingly. You can use the same scaling factor L to scale Ds, i.e. Ds2 = Ds(x2) = Ds(x/L).

3. Scale the reaction parameters (Vmax and Km): The reaction parameters, Vmax and Km, are not dependent on the space domain and thus do not need to be scaled. They can remain the same as before.

4. Scale the initial conditions: As the initial conditions are specified in terms of x, they will also need to be scaled using the same factor L. So, the initial condition for S should be S(1,:) = cos(x2*pi/(2*Length)).

5. Scale the boundary conditions: Since you have a Neumann boundary condition at x = 0, it will not be affected by the scaling. However, for the boundary condition at x = 1, you will need to use the scaled value of x, i.e. S(j+1,numx) = S(j,numx)+dt*Ds*2*(S(j,numx-1)-S(j,numx))./dx.^2-((Vmax*dt*S(j,numx))/(Km+S(j,numx))).

6. Scale the exact solution: The exact solution should also be scaled using the same factor L, i.e. S_exact = cos(x2*pi/(2*Length))*exp(-(pi*sqrt(Ds2)/(2*Length))^2*tmax).

By following these steps, your code should be able to accommodate the scaled space domain. However, it is always recommended to test your code for different values of L to ensure its accuracy.