Euler and FD method

7 次查看（过去 30 天）
Rana Önem
Rana Önem 2021-12-30
评论： Rana Önem 2022-1-2
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
Torsten 2021-12-30
Why don't you use MATLAB's "pdepe" for a solution ?
Rana Önem
Rana Önem 2021-12-30
I need to use finite diffrence method

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

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

类别

在 Help Center 和 File Exchange 中查找有关 Calculus 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by