Asked by Noemi ZM
on 14 Jun 2019

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 by John D'Errico
on 15 Jun 2019

Edited by John D'Errico
on 15 Jun 2019

Accepted 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.

Noemi ZM
on 15 Jun 2019

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 ZM
on 15 Jun 2019

Sign in to comment.

Opportunities for recent engineering grads.

Apply Today
## 0 Comments

Sign in to comment.