Euler and FD method

Hi, I have to solve a diffrential equation and I write a code like this but at the line 37 I get an error "Subscript indices must either be real positive integers or logicals." How can I solve this? (Line 37 is written in bold)

clear all
clc

% p=neutron flux function, macroscopic cross section of fission*nu=b,
% macroscopic crossection of absorption=c 
% 1/v* dp(x,t)/dt=b*p(x,t)- c*p(x,t)-D*d^2p(x,t)/dx^2)+ S(x,t) 
% p= neutron flux= p(x,t) 
% dp(x,t)/dt (partial derivatives) = p(x,t+dt)-p(x,t)/dt first derivative
% d^2p(x,t)/dx^2) (second derivatives)= p(x-dx,t)-2*p(x,t)+p(x+dx,t)/dx^2
%%% 1/v*(p(x,t+dt)-p(x,t)/dt)=b*p(x,t)- c*p(x,t)-D*p(x-dx,t)-2*p(x,t)+p(x+dx,t)/dx^2
% Solve for p(x,t+dt) 
%%% p(x,t+dt)=v*(b*p(x,t)- c*p(x,t)-D*p(x-dx,t)-2*p(x,t)+p(x+dx,t)/dx^2)+p(x,t)/dt

v=10^4; % cm/s 
D=2; %cm 
c=0.099; %cm^-1 
b=0.097; %cm^-1 
dx=1; %cm Diffrence for x 
dt=0.248; 
a=50; %cm 
x=-a:a;
N=2*a/dx; %number of elements
M=200; %number of elements 
xvec=1:dx:N-1;
tvec=1:dt:M;
%What is p function?
pmat=zeros(length(xvec),length(tvec));
%İnitial condition
t0=1;
pmat(xvec,t0)=0;
%Boundary condition
pmat(a,:)=0;
S0=10;
S(xvec)=S0*cos(pi*xvec/2*a);
for tdt=1:length(tvec)-1
    for idx=1:length(xvec)-1
        pmat(idx,tdt+1)= v*(b*pmat(idx,tdt)- c*pmat(idx,tdt)-D*pmat(idx-1,tdt)-2*pmat(idx,tdt)+pmat(idx+1,tdt)/dx^2)+pmat(idx,tvec)/tdt+(S0*cos(pi*idx/2*a));
    end 
end 

[Torsten]

Why don't you use MATLAB's "pdepe" for a solution ?

[Rana Önem]

I need to use finite diffrence method

Answer 1

[Torsten]

function main
v=10^4; % cm/s 
D=2; %cm
c=0.099; %cm^-1 
b=0.097; %cm^-1 
dx=1; %cm Diffrence for x 
dt=0.248; 
a=50; %cm 
x=-a:dx:a;
N=2*a/dx+1; %number of elements
M=200; %number of elements 
%What is p function?
pmat=zeros(N,M);
%İnitial condition
pmat(:,1)=0;
%Boundary condition
pmat(1,:)=0;   % inserted to make the code work
% a boundary condition for p at x=-a is missing
pmat(N,:)=0;
S0=10;
S=S0*cos(pi*x/(2*a));
for tdt=1:M-1
    for idx=2:N-1
        pmat(idx,tdt+1)= pmat(idx,tdt) + v*dt*(b*pmat(idx,tdt)- c*pmat(idx,tdt)-D*(pmat(idx-1,tdt)-2*pmat(idx,tdt)+pmat(idx+1,tdt))/dx^2+S(idx));
    end 
end 
plot(x.',pmat(:,10))
end

[Rana Önem]

Thank you so much!

Answer 2

[John D'Errico]

What is dt?

dt=0.248; 

What is tvec? Is tvec composed of integers? (No, because you used the colon operator with a stride of 0.248.).

tvec=1:dt:M;

What are you then doing with tvec? You index into pmat. That is, in ONE place in this line:

pmat(idx,tdt+1)= v*(pmat(idx,tdt) + dt*(b*pmat(idx,tdt)- c*pmat(idx,tdt)-D*(pmat(idx-1,tdt)-2*pmat(idx,tdt)+pmat(idx+1,tdt))/dx^2+S(idx)));

you use tvec as an index into pmat.

pmat(idx,tvec)

I might giuess you did not want to do that. But how can I know? ;-)

When you get an error message, read what it tells you. Think about the possibility that you are using an index that is not an integer value.