In-Situ Combustion Model - MATLAB code

I need to solve in MATLAB the following system:

[![][1]][1]

And for it I would like to use these datas and scheme (Crank-Nicolson)

With this initial condition:

I really do not know how I can do this on MATLAB. I tried to put on MATLAB all my datas, as you can see but I do not know as I can finish.

About this code I have got the following error. This error arise in all the pharenteses of the scheme:

Follow the code. I think I need to add some informartion about approximations, for instance, for

. But I don't know exactly where to put these and it seems that these are not about the error that I get.

Many thanks for any help!

tic
    clear all
    close all
    %% INPUT PARAMETERS
    Xmax=0.05; % x maximal value
    Xmin=0; % x maximal value
    tmax=0.008; % x maximal value
    tmin=0; % x maximal value
    N=101 % number of iterations
    dt=10^(-5); % time step value
    tol=10^(-6); % precison on the staionnary solution
    Pet=1406;
    beta=7.44*10^(10);
    epsilon=93.8;
    theta0=3.67;
    u=3.76;
    phi=beta*(1-eta)*e^(-epsilon/(theta+theta0));
    rho=theta0/(theta+theta0);
    V0=[theta;eta];
    % discretise the domain
    dx= (xmax-xmin)/N;

    %%
    % Initial conditions
    theta(0,t)=0;
    eta(0,t)=1;
    theta(x,0)=0;
    eta(x,0)=0;
    
    %%
    V=V0;
    VCN=V
    
    %% 
    % loop through time
    nsteps= tmax/dt;
     for n=1 : nsteps
        % calculate the scheme
        for i = 1 : N+1
            for j=1 : N+1
            
            theta(x(j),t(i+1))+u*dt/(4*dx)(eta(x(j+1),t(i+1))*theta(x(j+1),t(i+1))-eta(x(j-1),t(i+1))*theta(x(j-1),t(i+1)))- 
            dt/(2*Pet*dx)*(theta(x(j-1),t(i+1))-theta(x(j),t(i+1))+theta(x(j+1),t(i+1)))-dt/2*phi(x(j),t(i+1))==theta(x(j),t(i))-
            u*dt/(4*dx)*(eta(x(j+1),t(i))*theta(x(j+1),t(i))-eta(x(j-1),t(i))*theta(x(j-1),t(i)))+dt/(Pet*2dx)*(theta(x(j-1),t(i))-
            2*theta(x(j),t(i))+theta(x(j+1),t(i)))+dt/2*phi(x(j),t(i));
            
            dt/2*phi(x(j),t(i+1))+eta(x(j),t(i+1))==eta(x(j),t(i))-dt/2*phi(x(j),t(i));
            end
        end
    
        % update t and V
        t=t+dt;
        V=VCN;
     end 
     
    function x=x(j);
    x(j)=j*dx 
    end
    
    function t=t(i);
    t(i)=i*dt 
    end
    
    function eta=eta(x,t);
    end
    
    function theta=theta(x,t);
    end
    
    function phi=phi(x,t);
    phi(x,t)=beta*(1-eta(x,t))*e^(-epsilon/(theta(x,t)+theta0))
    end 
    
        % update t and V
        t=t+dt;
        V=VCN;
    
    % plot solution
        set(fig,'YData',V(1,:));
        set(tit,'String',strcat('n =',num2str(n)));
        drawnow;
        
     end        

Answer 1

[John D'Errico]

This is a simple error of syntax, which seems a bit strange for someone who has written such a long code to misunderstand. You should spend a little time learning the basic syntax of MATLAB.

MATLAB uses * to multiply. Or .* if you want to make sure that this is an element-wise matiplication. But it uses an operator to multiply.

The point is though, that

(x-y)(x+y)

is NOT the product of those two numbers, at least not in MATLAB. It may be that when you write it on paper, etc. But in MATLAB? What I wrote there is a syntax error. If you want to multiply numbers, then you genrally need to write

(x-y).*(x+y)

in case they are vectors or matrices. If both x and y are scalars, then it would be sufficient to write:

(x-y)*(x+y)

But in any case, you can see that I put an operator in there, between the pieces. Now go back and read the error message. Look at the point that is indicated. MATLAB sees a situation just like what I wrote, and it does not understand what you intended there, because what you wrote is not proper MATLAB syntax. For example, when MATLAB sees this generic form:

A(B)

it assumes that A is a function, that will be called with an operand B. So, when I write

(A)(B)

it gets confused. Is A a function, and B an argument? Computers need syntax to communicate with you.

[Noemi Zeraick Monteiro]

    tic
    clear all
    close all
    %% INPUT PARAMETERS
    xmax=0.05; % x maximal value
    xmin=0; % x maximal value
    tmax=0.008; % x maximal value
    tmin=0; % x maximal value
    dt=10^(-5); % time step value
    tol=10^(-6); % precison on the staionnary solution
    Pet=1406;
    beta=7.44*10^(10);
    epsilon=93.8;
    theta0=3.67;
    u=3.76;
    phi==beta*(1-eta)*e^(-epsilon/(theta+theta0));
    rho==theta0/(theta+theta0);
    V0==[theta;eta];
    % discretise the domain
    dX=(xmax-xmin)/500;
    
    
    %%
    V=V0;
    VCN=V
    %% 
    
    % Initial conditions
    for i=0: 800;
        if x==0; 
        theta(x,t(i))==0;
        eta(x,t(i))==1;
        end 
    end
    
    for i=0: 800;
        if t==0;
        theta(x(i),t)==0;
        eta(x(i),t)==0;
        end 
    end
    
    % loop through time
    nsteps= tmax/dt;
     for n=1 : nsteps
    
    
        % calculate the scheme
        for i = 1 : 800
            for j=1 : 800
            
               theta(x(j),t(i+1))+u*dt*(1/(4*dX))*(rho(x(j+1),t(i+1))*theta(x(j+1),t(i+1))-rho(x(j-1),t(i+1))*theta(x(j-1),t(i+1)))-(dt/(2*Pet*dX*dX))*(theta(x(j-1),t(i+1))-theta(x(j),t(i+1))+theta(x(j+1),t(i+1)))-(dt/2)*phi(x(j),t(i+1))== theta(x(j),t(i))-u*dt*(1/(4*dX))*(rho(x(j+1),t(i))*theta(x(j+1),t(i))-rho(x(j-1),t(i))*theta(x(j-1),t(i)))+ (dt/(Pet*2*dX*dX))*(theta(x(j-1),t(i))-2*theta(x(j),t(i))+theta(x(j+1),t(i)))+(dt/2)*phi(x(j),t(i));
    
    
               (dt/2)*phi(x(j),t(i+1))+eta(x(j),t(i+1))== eta(x(j),t(i))-(dt/2)*phi(x(j),t(i));
            
            end
        end
    
        % update V
        
        V=VCN;
     end 
     
     
     % plot solution
        set(fig,'YData',V(1,:));
        set(tit,'String',strcat('n =',num2str(n)));
        drawnow;
       
    
        
    function x=x(j);
    x(j)=j*dx 
    end
    
    function t=t(i);
    t(i)=i*dt 
    end
    
    function eta=eta(x,t);
    end
    
    function theta=theta(x,t);
    end
    
    function phi=phi(x,t);
    phi(x,t)=beta*(1-eta(x,t))*e^(-epsilon/(theta(x,t)+theta0));
    end 

[Noemi Zeraick Monteiro]

Hi. Sorry, I've got a hard activity and I don't have experience on MATLAB, indeed my deadline is hard, so I'm sorry for my mistakes. I've tried again I got the following error.