MATLAB: In-Situ Combustion Model – MATLAB code

mathematicsnumerical method

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

Best Answer

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.
Related Question