Euler and FD method

7 Ansichten (letzte 30 Tage)
Rana Önem
Rana Önem am 30 Dez. 2021
Kommentiert: Rana Önem am 2 Jan. 2022
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
  2 Kommentare
Torsten
Torsten am 30 Dez. 2021
Why don't you use MATLAB's "pdepe" for a solution ?
Rana Önem
Rana Önem am 30 Dez. 2021
I need to use finite diffrence method

Melden Sie sich an, um zu kommentieren.

Akzeptierte Antwort

Torsten
Torsten am 30 Dez. 2021
Bearbeitet: Torsten am 30 Dez. 2021
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

Weitere Antworten (1)

John D'Errico
John D'Errico am 30 Dez. 2021
Bearbeitet: John D'Errico am 30 Dez. 2021
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.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by