Euler and FD method
조회 수: 7 (최근 30일)
이전 댓글 표시
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
채택된 답변
Torsten
2021년 12월 30일
편집: Torsten
2021년 12월 30일
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
추가 답변 (1개)
John D'Errico
2021년 12월 30일
편집: John D'Errico
2021년 12월 30일
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.
댓글 수: 0
참고 항목
카테고리
Help Center 및 File Exchange에서 Calculus에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!