Why do I have an error in pdepe solving?

1 次查看（过去 30 天）
Rahul
Rahul 2023-10-18
编辑： Torsten 2023-10-18
Ran in:
Hi all,
I'm having an error while solving a pdepe. It gives me the error
"Index exceeds the number of array elements. Index must not exceed 4."
It will help me immensely if someone kindly help me in understanding as why is it here.
I know in general it happens if array exceeds the no.
But precisely, this seems not the case here in mine.
kindly find the code as attached.
pde2fshear_v4Perturbed_randomWalk()
drift_vel = -0.1000 dCoeff = 0.0300
Index exceeds the number of array elements. Index must not exceed 4.

Error in solution>pdex2pde (line 544)
s(4) = chi_growth*(term1*theta_heaviside1-gamma_nonlin*u(5)^alpha_nonlin-alpha_0*u(3)-alpha_v*u(4))*u(5) -drift_vel*DuDx(5);

Error in pdepe (line 242)
[c,f,s] = feval(pde,xi(1),t(1),U,Ux,varargin{:});

Error in solution>pde2fshear_v4Perturbed_randomWalk (line 152)
sol = pdepe(m,@pdex2pde,@pdex2ic,@pdex2bc,x,t);
% ****************************************************************
% * Program: pde2fshear_v4Perturbed (Random walk model) *
% * Description: Solving coupled thermal, particle transport & *
% * turbulence intensity equations, suppressed by flow shear
% * and zonal shear *
% * Author: Rahul Chakraborty *
% * Email: 6310230015@email.psu.ac.th *
% * Date: 2023-04-29 *
% * Version: 1.0 *
% ****************************************************************
% solve 3-F bifurcation model
function pde2fshear_v4Perturbed_randomWalk
global chi0; % declare global variables
global D0;
global chi1;
global D1;
global alpha_chi;
global alpha_D;
global H0;
global S0;
global H;
global S;
global data;
global track;
global xstep;
global tstep;
global count;
global theta_heaviside1;
global chi_growth;
global beta_growth; % strength of non-linear growth
global delta_suppress;
global lambda_suppress;
global v_e;
global chi_ano;
global D_ano;
global sigma_turb;
%constants for random walk model
global gamma_nonlin;
global alpha_nonlin;
global drift_vel;
global dCoeff;
global rho_0;
global gamma_damp;
global alpha_0;
global alpha_v;
%define values for constants
data.constant.critgradpressure = 1.2;
data.constant.critgraddensity = 1.0;
count = 1;
chi0 = 0.5; D0 = 0.5; chi1 = 5.0; D1 = 5.0;
alpha_chi = 0.1; alpha_D = 0.1;
H0 = 25; S0 = 21;
track = 1;
chi_growth=0.5;
lambda_suppress=0.01;
sigma_turb=1.5;
gamma_nonlin = 10;
alpha_nonlin = 1;
rho_0=1;
gamma_damp=5;
alpha_0=5;
alpha_v=1;
%position and time grids information
xstep = 100;
tstep = 100;
xmin = 0;
xmax = 1;
tmin = 0;
tmax = 10;
m = 0; %define type of equation to solve
%Preallocate vectors for speed improvement
grad_u1 = zeros(tstep,xstep);
grad_u2 = zeros(tstep,xstep);
grad_u3 = zeros(tstep,xstep);
grad_u4 = zeros(tstep,xstep);
grad_u5 = zeros(tstep,xstep);
curve_u1 = zeros(tstep,xstep);
curve_u2 = zeros(tstep,xstep);
curve_u3 = zeros(tstep,xstep);
curve_u4 = zeros(tstep,xstep);
curve_u5 = zeros(tstep,xstep);
flowshear_p = zeros(tstep,xstep);
flowshear_n = zeros(tstep,xstep);
Q = zeros(tstep,xstep);
Q0 = zeros(tstep,xstep);
Q1 = zeros(tstep,xstep);
neo_p = zeros(tstep,xstep);
ano_p = zeros(tstep,xstep);
Gam = zeros(tstep,xstep);
Gam0 = zeros(tstep,xstep);
Gam1 = zeros(tstep,xstep);
neo_n = zeros(tstep,xstep);
ano_n = zeros(tstep,xstep);
%define first inner half of the plasma
%x = linspace(xmin,xmax/2,xstep/5);
x = linspace(xmin,xmax,xstep);
t = linspace(tmin,tmax,tstep);
%Add smaller grid size near plasma edge
%for i=(xstep/5)+1:xstep
% x(i) = x(i-1)+(xmax/2)/(xstep-(xstep/5));
%end
data.variable.x = x;
%Turbulence intensity equation for Random Walk
numberofWalk = 5;
stepsperWalk = 5;
std_step = 0.1;
mean_wait = 1.0; %set as appropriate
delta_X = sign(randn(1,stepsperWalk)) * std_step;
delta_T = sign(rand(1,stepsperWalk)) * mean_wait;
t_total = transpose(cumsum([0, delta_T]));
x_total = transpose(cumsum([0, delta_X]));
stem(t_total, x_total)
xlabel('time'); ylabel('displacement');
xt = zeros(stepsperWalk,2);
for walk = 1:numberofWalk
for step = 1:stepsperWalk
xt(step+1,1) = xt(step,1) + x_total(step);
xt(step+1,2) = xt(step,2) + t_total(step);
end
end
plot(xt(:,2),xt(:,1))
drift_vel = mean(xt(step,1)/xt(step,2));
dCoeff = mean(xt(step,1)^2/(2*xt(step,2)));
disp('drift_vel = '); disp(drift_vel);
disp('dCoeff = '); disp(dCoeff);
sol = pdepe(m,@pdex2pde,@pdex2ic,@pdex2bc,x,t);
% Extract the first solution component as u1 = pressure
% second solution component as u2 = density
% third solution component as u3 = zonal flow
% fourth as u4 = mean flow
% fifth as u5 = turb intensity
u1 = sol(:,:,1);
u2 = sol(:,:,2);
u3 = sol(:,:,3);
u4 = sol(:,:,4);
u5 = sol(:,:,5);
%grad_u = gradient(u,(x(2)-x(1)));
for j=1:tstep
for i = 1:xstep/5
if i == 1
grad_u1(j,i) = (u1(j,2)-u1(j,1))/(x(2)-x(1));
grad_u2(j,i) = (u2(j,2)-u2(j,1))/(x(2)-x(1));
grad_u3(j,i) = (u3(j,2)-u3(j,1))/(x(2)-x(1));
grad_u4(j,i) = (u4(j,2)-u4(j,1))/(x(2)-x(1));
grad_u5(j,i) = (u5(j,2)-u5(j,1))/(x(2)-x(1));
elseif i == xstep/5
grad_u1(j,i) = (u1(j,i)-u1(j,i-1))/(x(i)-x(i-1));
grad_u2(j,i) = (u2(j,i)-u2(j,i-1))/(x(i)-x(i-1));
grad_u3(j,i) = (u3(j,i)-u3(j,i-1))/(x(i)-x(i-1));
grad_u4(j,i) = (u4(j,i)-u4(j,i-1))/(x(i)-x(i-1));
grad_u5(j,i) = (u5(j,i)-u5(j,i-1))/(x(i)-x(i-1));
else
grad_u1(j,i) = (u1(j,i+1)-u1(j,i-1))/(x(i+1)-x(i-1));
grad_u2(j,i) = (u2(j,i+1)-u2(j,i-1))/(x(i+1)-x(i-1));
grad_u3(j,i) = (u3(j,i+1)-u3(j,i-1))/(x(i+1)-x(i-1));
grad_u4(j,i) = (u4(j,i+1)-u4(j,i-1))/(x(i+1)-x(i-1));
grad_u5(j,i) = (u5(j,i+1)-u5(j,i-1))/(x(i+1)-x(i-1));
end
end
for i=(xstep/5)+1:xstep
if i == xstep/5+1
grad_u1(j,i) = (u1(j,i+1)-u1(j,i))/(x(i+1)-x(i));
grad_u2(j,i) = (u2(j,i+1)-u2(j,i))/(x(i+1)-x(i));
grad_u3(j,i) = (u3(j,i+1)-u3(j,i))/(x(i+1)-x(i));
grad_u4(j,i) = (u4(j,i+1)-u4(j,i))/(x(i+1)-x(i));
grad_u5(j,i) = (u5(j,i+1)-u5(j,i))/(x(i+1)-x(i));
elseif i == xstep
grad_u1(j,i) = (u1(j,i)-u1(j,i-1))/(x(i)-x(i-1));
grad_u2(j,i) = (u2(j,i)-u2(j,i-1))/(x(i)-x(i-1));
grad_u3(j,i) = (u3(j,i)-u3(j,i-1))/(x(i)-x(i-1));
grad_u4(j,i) = (u4(j,i)-u4(j,i-1))/(x(i)-x(i-1));
grad_u5(j,i) = (u5(j,i)-u5(j,i-1))/(x(i)-x(i-1));
else
grad_u1(j,i) = (u1(j,i+1)-u1(j,i-1))/(x(i+1)-x(i-1));
grad_u2(j,i) = (u2(j,i+1)-u2(j,i-1))/(x(i+1)-x(i-1));
grad_u3(j,i) = (u3(j,i+1)-u3(j,i-1))/(x(i+1)-x(i-1));
grad_u4(j,i) = (u4(j,i+1)-u4(j,i-1))/(x(i+1)-x(i-1));
grad_u5(j,i) = (u5(j,i+1)-u5(j,i-1))/(x(i+1)-x(i-1));
end
end
end
for i=1:tstep
for j=1:xstep
v_e = -grad_u1(i,j)*grad_u2(i,j)/u2(i,j)^2; % -g_p*g_n/n^2
flowshear_p(i,j) = 1+ alpha_chi*v_e^2;
flowshear_n(i,j) = 1+ alpha_D*v_e^2;
if abs(grad_u1(i,j)) < abs(data.constant.critgradpressure) %&& abs(grad_u2(i,j)) < abs(data.constant.critgraddensity)
Q(i,j) = -grad_u1(i,j)*chi0;
Q0(i,j) = Q(i,j);
Q1(i,j) = 0;
neo_p(i,j) = chi0*(1+grad_u1(i,j))/(1+grad_u1(i,j));
ano_p(i,j) = 0;
Gam(i,j) = -grad_u2(i,j)*D0;
Gam0(i,j) = Gam(i,j);
Gam1(i,j) = 0;
neo_n(i,j) = D0*(1+grad_u2(i,j))/(1+grad_u2(i,j));
ano_n(i,j) = 0;
elseif abs(grad_u1(i,j)) >= abs(data.constant.critgradpressure) %&& abs(grad_u2(i,j)) < abs(data.constant.critgraddensity)
Q(i,j) = (chi0*(-grad_u1(i,j)) + chi1*(-grad_u1(i,j)+data.constant.critgradpressure)*u5(i,j)/flowshear_p(i,j))*(-grad_u1(i,j));
Q0(i,j) = -grad_u1(i,j)*chi0;
Q1(i,j) = (chi1*(-grad_u1(i,j)+data.constant.critgradpressure)*u5(i,j)/flowshear_p(i,j))*(-grad_u1(i,j));
neo_p(i,j) = chi0*(1+grad_u1(i,j))/(1+grad_u1(i,j));
ano_p(i,j) = chi1*(abs(grad_u1(i,j))+data.constant.critgradpressure)*u5(i,j)/flowshear_p(i,j)*(-grad_u1(i,j));
Gam(i,j) = -grad_u2(i,j)*D0;
Gam0(i,j) = Gam(i,j);
Gam1(i,j) = 0;
neo_n(i,j) = D0*(1+grad_u2(i,j))/(1+grad_u2(i,j));
ano_n(i,j) = 0;
else
Q(i,j) = (chi0*(-grad_u1(i,j)) + chi1*(-grad_u1(i,j)+data.constant.critgradpressure)*u5(i,j)/flowshear_p(i,j))*(-grad_u1(i,j));
Q0(i,j) = -grad_u1(i,j)*chi0;
Q1(i,j) = (chi1*(-grad_u1(i,j)+data.constant.critgradpressure)*u5(i,j)/flowshear_p(i,j))*(-grad_u1(i,j));
neo_p(i,j) = chi0*(1+grad_u1(i,j))/(1+grad_u1(i,j));
ano_p(i,j) = chi1*(abs(grad_u1(i,j))+data.constant.critgradpressure)*u5(i,j)/flowshear_p(i,j)*(-grad_u1(i,j));
Gam(i,j) = -grad_u2(i,j)*(D0 + D1*(-grad_u2(i,j)+data.constant.critgraddensity)*u5(i,j)/flowshear_n(i,j));
Gam0(i,j) = -grad_u2(i,j)*D0;
Gam1(i,j) = -grad_u2(i,j)*D1*(-grad_u2(i,j)+data.constant.critgraddensity)*u5(i,j)/flowshear_n(i,j);
neo_n(i,j) = D0*(1+grad_u2(i,j))/(1+grad_u2(i,j));
ano_n(i,j) = D1*(abs(grad_u2(i,j))+data.constant.critgraddensity)/flowshear_n(i,j);
end
end
end
%curve_u = gradient(grad_u,(x(2)-x(1)));
for j=1:tstep
for i = 1:xstep/5
if i == 1
curve_u1(j,i) = (grad_u1(j,2)-grad_u1(j,1))/(x(2)-x(1));
curve_u2(j,i) = (grad_u2(j,2)-grad_u2(j,1))/(x(2)-x(1));
curve_u3(j,i) = (grad_u3(j,2)-grad_u3(j,1))/(x(2)-x(1));
curve_u4(j,i) = (grad_u4(j,2)-grad_u4(j,1))/(x(2)-x(1));
curve_u5(j,i) = (grad_u5(j,2)-grad_u5(j,1))/(x(2)-x(1));
elseif i == xstep/5
curve_u1(j,i) = (grad_u1(j,i)-grad_u1(j,i-1))/(x(i)-x(i-1));
curve_u2(j,i) = (grad_u2(j,i)-grad_u2(j,i-1))/(x(i)-x(i-1));
curve_u3(j,i) = (grad_u3(j,i)-grad_u3(j,i-1))/(x(i)-x(i-1));
curve_u4(j,i) = (grad_u4(j,i)-grad_u4(j,i-1))/(x(i)-x(i-1));
curve_u5(j,i) = (grad_u5(j,i)-grad_u5(j,i-1))/(x(i)-x(i-1));
else
curve_u1(j,i) = (grad_u1(j,i+1)-grad_u1(j,i-1))/(x(i+1)-x(i-1));
curve_u2(j,i) = (grad_u2(j,i+1)-grad_u2(j,i-1))/(x(i+1)-x(i-1));
curve_u3(j,i) = (grad_u3(j,i+1)-grad_u3(j,i-1))/(x(i+1)-x(i-1));
curve_u4(j,i) = (grad_u4(j,i+1)-grad_u4(j,i-1))/(x(i+1)-x(i-1));
curve_u5(j,i) = (grad_u5(j,i+1)-grad_u5(j,i-1))/(x(i+1)-x(i-1));
end
end
for i=(xstep/5)+1:xstep
if i == xstep/5+1
curve_u1(j,i) = (grad_u1(j,i+1)-grad_u1(j,i))/(x(i+1)-x(i));
curve_u2(j,i) = (grad_u2(j,i+1)-grad_u2(j,i))/(x(i+1)-x(i));
curve_u3(j,i) = (grad_u3(j,i+1)-grad_u3(j,i))/(x(i+1)-x(i));
curve_u4(j,i) = (grad_u4(j,i+1)-grad_u4(j,i))/(x(i+1)-x(i));
curve_u5(j,i) = (grad_u5(j,i+1)-grad_u5(j,i))/(x(i+1)-x(i));
elseif i == xstep
curve_u1(j,i) = (grad_u1(j,i)-grad_u1(j,i-1))/(x(i)-x(i-1));
curve_u2(j,i) = (grad_u2(j,i)-grad_u2(j,i-1))/(x(i)-x(i-1));
curve_u3(j,i) = (grad_u3(j,i)-grad_u3(j,i-1))/(x(i)-x(i-1));
curve_u4(j,i) = (grad_u4(j,i)-grad_u4(j,i-1))/(x(i)-x(i-1));
curve_u5(j,i) = (grad_u5(j,i)-grad_u5(j,i-1))/(x(i)-x(i-1));
else
curve_u1(j,i) = (grad_u1(j,i+1)-grad_u1(j,i-1))/(x(i+1)-x(i-1));
curve_u2(j,i) = (grad_u2(j,i+1)-grad_u2(j,i-1))/(x(i+1)-x(i-1));
curve_u3(j,i) = (grad_u3(j,i+1)-grad_u3(j,i-1))/(x(i+1)-x(i-1));
curve_u4(j,i) = (grad_u4(j,i+1)-grad_u4(j,i-1))/(x(i+1)-x(i-1));
curve_u5(j,i) = (grad_u5(j,i+1)-grad_u5(j,i-1))/(x(i+1)-x(i-1));
end
end
end
%to save parameters and variables
data.constant.chi0 = chi0;
data.constant.D0 = D0;
data.constant.chi1 = chi1;
data.constant.D1 = D1;
data.constant.alphachi = alpha_chi;
data.constant.alphaD = alpha_D;
data.constant.H0 = H0;
data.constant.S0 = S0;
data.variable.pressure = u1;
data.variable.density = u2;
data.variable.zonalflow = u3;
data.variable.meanflow = u4;
data.variable.turbulenceintensity = u5;
data.variable.gradpressure = grad_u1;
data.variable.graddensity = grad_u2;
data.variable.gradzonalflow = grad_u3;
data.variable.gradmeanflow = grad_u4;
data.variable.gradturbulenceintensity = grad_u5;
data.variable.curvepressure = curve_u1;
data.variable.curvedensity = curve_u2;
data.variable.curvezonalflow = curve_u3;
data.variable.curvemeanflow = curve_u4;
data.variable.curveturbulenceintensity = curve_u5;
data.variable.x = x;
data.variable.t = t;
data.variable.Q = Q;
data.variable.Gamma = Gam;
data.variable.Q0 = Q0;
data.variable.Gamma0 = Gam0;
data.variable.neo_P = neo_p;
data.variable.neo_n = neo_n;
data.variable.Q1 = Q1;
data.variable.Gam1 = Gam1;
data.variable.ano_p = ano_p;
data.variable.ano_n = ano_n;
data.variable.heatsource = H;
data.variable.particlesource = S;
data.variable.wexb_p = flowshear_p;
data.variable.wexb_n = flowshear_n;
%data.predicted.Q = Q_pre;
%data.predicted.Gam = Gam_pre;
%data.predicted.gradientp = pp_pre;
%data.predicted.gradientn = nn_pre;
data.control.xgrid = xstep;
data.control.tgrid = tstep;
data.control.xmin = xmin;
data.control.xmax = xmax;
data.control.tmin = tmin;
data.control.tmax = tmax;
figure;
subplot(2,1,1,'FontSize',22)
plot(x,u1(end,:),'.');
axis ([0 1 0 20]);
ylabel('p')
xlabel('r/a')
%set(gca,'XTick',[0,0.2,0.4,0.6,0.8])
head = strcat('Time = ',num2str(t(end)),' s');
title(head);
grid on
subplot(2,1,2,'FontSize',22)
plot(x,grad_u1(end,:),'.');
axis([0 1 -60 0]);
ylabel('p\prime')
xlabel('r/a')
grid on
figure;
subplot(1,2,1,'FontSize',22)
plot(x,curve_u1(end,:),'.');
axis ([0 1 -1400 0]);
ylabel('p\prime\prime')
xlabel('r/a')
%set(gca,'XTick',[0,0.2,0.4,0.6,0.8])
head = strcat('Time = ',num2str(t(end)),' s');
title(head);
grid on
subplot(1,2,2,'FontSize',22)
plot(abs(grad_u1(end,:)),Q(end,:),'.');
%plot(pp_pre,Q_pre,'color','blue','LineWidth',3);
%hold on
%axis ([0 60 0 60]);
ylabel('Q')
xlabel('|p\prime|')
%title('Bifurcation diagram')
grid on
figure;
subplot(2,1,1,'FontSize',22)
plot(x,u2(end,:),'.');
axis ([0 1 0 20]);
ylabel('n')
xlabel('r/a')
%set(gca,'XTick',[0,0.2,0.4,0.6,0.8])
head = strcat('Time = ',num2str(t(end)),' s');
title(head);
grid on
subplot(2,1,2,'FontSize',22)
plot(x,grad_u2(end,:),'.');
axis([0 1 -60 0]);
ylabel('n\prime')
xlabel('r/a')
grid on
figure;
subplot(1,2,1,'FontSize',22)
plot(x,curve_u2(end,:),'.');
axis ([0 1 -1400 0]);
ylabel('n\prime\prime')
xlabel('r/a')
%set(gca,'XTick',[0,0.2,0.4,0.6,0.8])
head = strcat('Time = ',num2str(t(end)),' s');
title(head);
grid on
subplot(1,2,2,'FontSize',22)
plot(abs(grad_u2(end,:)),Gam(end,:),'.');
%plot(nn_pre,Gam_pre,'color','blue','LineWidth',3);
%hold on
%axis ([0 60 0 60]);
ylabel('\Gamma')
xlabel('|n\prime|')
%title('Bifurcation diagram')
grid on
figure;
subplot(2,1,1,'FontSize',22)
plot(x,u3(end,:),'.');
axis ([0 1 0 20]);
ylabel('I')
xlabel('r/a')
%set(gca,'XTick',[0,0.2,0.4,0.6,0.8])
head = strcat('Time = ',num2str(t(end)),' s');
title(head);
grid on
subplot(2,1,2,'FontSize',22)
plot(x,grad_u3(end,:),'.');
axis([0 1 -60 0]);
ylabel('I\prime')
xlabel('r/a')
grid on
figure;
subplot(1,2,1,'FontSize',22)
plot(x,curve_u3(end,:),'.');
axis ([0 1 -1400 0]);
ylabel('I\prime\prime')
xlabel('r/a')
%set(gca,'XTick',[0,0.2,0.4,0.6,0.8])
head = strcat('Time = ',num2str(t(end)),' s');
title(head);
grid on
subplot(1,2,2,'FontSize',22)
plot(abs(grad_u3(end,:)),Q(end,:),'.');
%plot(pp_pre,Q_pre,'color','blue','LineWidth',3);
%hold on
%axis ([0 60 0 60]);
ylabel('Q')
xlabel('|I\prime|')
%title('Bifurcation diagram')
grid on
end
% --------------------------------------------------------------
function [c,f,s] = pdex2pde(x,t,u,DuDx)
global chi0;
global D0;
global chi1;
global D1;
global H0;
global S0;
global data;
global track;
global xstep;
global alpha_chi;
global alpha_D;
global chi_growth; % total growth rate
global beta_growth; % strength of non-linear growth
global delta_suppress;
global length;
global tau_conf;
global I0; % edge localized source of turbulence intensity
global theta_heaviside1;
%global count;
global v_e;
global drift_vel;
%constants for random walk model
global gamma_nonlin;
global alpha_nonlin;
global drift_vel;
global dCoeff;
global alpha_0;
global rho_0;
global gamma_damp;
global alpha_v;
c = [1;1;1;1];
length=1;
v_e = -DuDx(1)*DuDx(2)/u(2)^2; % -g_p*g_n/n^2
flowshear_p = 1+ alpha_chi*v_e^2;
flowshear_n = 1+ alpha_D*v_e^2;
u(4) = mean(v_e)^2;
term1 = abs(DuDx(1))-data.constant.critgradpressure;
% Implementing Heaviside function for H-mode in p, n and I equations
if term1 > 0
theta_heaviside1=1;
else
theta_heaviside1=0;
end
s(3) = -gamma_damp*u(3);
s(4) = chi_growth*(term1*theta_heaviside1-gamma_nonlin*u(5)^alpha_nonlin-alpha_0*u(3)-alpha_v*u(4))*u(5) -drift_vel*DuDx(5);
s = [H0*exp(-100*x^2/length)+H0/2; S0*exp(-100*(x-0.9)^2/length)+S0/2; s(3);s(4)];
f = [chi0+chi1*u(5)/flowshear_p; D0+D1*u(5)/flowshear_n; alpha_0*u(3)*u(5)/(1+rho_0*u(4)); dCoeff].*DuDx; % flux term for CTRW
end
% --------------------------------------------------------------
function u0 = pdex2ic(x)
%u0 = [eps; eps; eps];
%u0 = [0.01; 0.01; 0.1*exp(-100*(x-1)^2)];
u0 = [0.4*(1-x^2); 0.4*(1-x^2);0.4*(1-x^2);0.1*exp(-100*(x-1)^2)];
end
% --------------------------------------------------------------
function [pl,ql,pr,qr] = pdex2bc(xl,ul,xr,ur,t)
pl = [0; 0; 0; 0];
ql = [1; 1; 1; 1];
pr = [ur(1); ur(2); ur(3); 0];
qr = [0.01; 0.001; 1; 1;];
% qr = [0; 0.1; 1];
%qr = [0; eps; 1];
end
%---------------------------------------------------------------
with regards,
rc

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

回答（1 个）

Shivam
Shivam 2023-10-18
Hi Rahul,
From the information provided, I understand that while executing the attached code, you are getting "index out of bound" issue.
You can debug in your code to know that the error is in line 544 because of indexing in vector 'u' for index 5 whereas the size of vector 'u' is 4x1 only.
>> u
u =
0.4000
0.4000
0.4000
0.0000
I hope this helps.
Thanks,
Shivam
  4 个评论
显示 2更早的评论
Rahul
Rahul 2023-10-18
编辑：Torsten 2023-10-18
Ran in:
Hi Torsten,
Thanks for your comment.
I did changes in the code. But still the error continues.
Attached herewith the code.Can u plz comment.
pde2fshear_v4Perturbed_randomWalk()
Warning: Failure at t=1.855077e-05. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (5.421011e-20) at time t.
Warning: Time integration has failed. Solution is available at requested time points up to t=0.000000e+00.
Index in position 1 exceeds array bounds. Index must not exceed 1.

Error in solution>pde2fshear_v4Perturbed_randomWalk (line 135)
grad_u1(j,i) = (u1(j,2)-u1(j,1))/(x(2)-x(1));
function pde2fshear_v4Perturbed_randomWalk
global chi0; % declare global variables
global D0;
global chi1;
global D1;
global alpha_chi;
global alpha_D;
global H0;
global S0;
global Y0;
global H;
global S;
global data;
global track;
global xstep;
global tstep;
global count;
global theta_heaviside1;
global chi_growth;
global lambda_suppress;
global v_e;
global chi_ano;
global D_ano;
global sigma_turb;
global gamma_nonlin;
global alpha_zero;
global alpha_v;
global rho_0;
global gamma_damp;
global alpha_nonlin;
global dCoeff;
%define values for constants
data.constant.critgradpressure = 1.2;
data.constant.critgraddensity = 1.0;
count = 1;
chi0 = 0.5;
D0 = 0.5;
chi1 = 0.5;
D1 = 0.5;
alpha_chi = 0.1;
alpha_D = 0.1;
H0 = 4;
S0 = 13;
track = 1;
chi_growth=20;
lambda_suppress=0.5 ;
sigma_turb=2.0;
chi_ano=10;
D_ano=10;
gamma_nonlin = 10;
gamma_damp = 1;
alpha_zero = 1;
alpha_v = 1;
rho_0 = 1;
alpha_nonlin = 1;
dCoeff = 1;
%position and time grids information
xstep = 100;
tstep = 100;
xmin = 0;
xmax = 1;
tmin = 0;
tmax = 10;
%tmax = 1;
m = 0; %define type of equation to solve
%Preallocate vectors for speed improvement
grad_u1 = zeros(tstep,xstep);
grad_u2 = zeros(tstep,xstep);
grad_u3 = zeros(tstep,xstep);
grad_u4 = zeros(tstep,xstep);
grad_u5 = zeros(tstep,xstep);
curve_u1 = zeros(tstep,xstep);
curve_u2 = zeros(tstep,xstep);
curve_u3 = zeros(tstep,xstep);
curve_u4 = zeros(tstep,xstep);
curve_u5 = zeros(tstep,xstep);
flowshear_p = zeros(tstep,xstep);
flowshear_n = zeros(tstep,xstep);
Q = zeros(tstep,xstep);
Q0 = zeros(tstep,xstep);
Q1 = zeros(tstep,xstep);
neo_p = zeros(tstep,xstep);
ano_p = zeros(tstep,xstep);
Gam = zeros(tstep,xstep);
Gam0 = zeros(tstep,xstep);
Gam1 = zeros(tstep,xstep);
neo_n = zeros(tstep,xstep);
ano_n = zeros(tstep,xstep);
%define first inner half of the plasma
%x = linspace(xmin,xmax/2,xstep/5);
x = linspace(xmin,xmax,xstep);
t = linspace(tmin,tmax,tstep);
%Add smaller grid size near plasma edge
%for i=(xstep/5)+1:xstep
% x(i) = x(i-1)+(xmax/2)/(xstep-(xstep/5));
%end
data.variable.x = x;
sol = pdepe(m,@pdex2pde,@pdex2ic,@pdex2bc,x,t);
% Extract the first solution component as u1 = pressure
% second solution component as u2 = density
% third solution component as u3 = turbulence intensity
% fourth solution component as u4 = zonal flow
% fifth solution component as u5 = mean flow
u1 = sol(:,:,1);
u2 = sol(:,:,2);
u3 = sol(:,:,3);
u4 = sol(:,:,4);
u5 = sol(:,:,5);
%grad_u = gradient(u,(x(2)-x(1)));
for j=1:tstep
for i = 1:xstep/5
if i == 1
grad_u1(j,i) = (u1(j,2)-u1(j,1))/(x(2)-x(1));
grad_u2(j,i) = (u2(j,2)-u2(j,1))/(x(2)-x(1));
grad_u3(j,i) = (u3(j,2)-u3(j,1))/(x(2)-x(1));
grad_u4(j,i) = (u4(j,2)-u4(j,1))/(x(2)-x(1));
grad_u5(j,i) = (u5(j,2)-u5(j,1))/(x(2)-x(1));
elseif i == xstep/5
grad_u1(j,i) = (u1(j,i)-u1(j,i-1))/(x(i)-x(i-1));
grad_u2(j,i) = (u2(j,i)-u2(j,i-1))/(x(i)-x(i-1));
grad_u3(j,i) = (u3(j,i)-u3(j,i-1))/(x(i)-x(i-1));
grad_u4(j,i) = (u4(j,i)-u4(j,i-1))/(x(i)-x(i-1));
grad_u5(j,i) = (u5(j,i)-u5(j,i-1))/(x(i)-x(i-1));
else
grad_u1(j,i) = (u1(j,i+1)-u1(j,i-1))/(x(i+1)-x(i-1));
grad_u2(j,i) = (u2(j,i+1)-u2(j,i-1))/(x(i+1)-x(i-1));
grad_u3(j,i) = (u3(j,i+1)-u3(j,i-1))/(x(i+1)-x(i-1));
grad_u4(j,i) = (u4(j,i+1)-u4(j,i-1))/(x(i+1)-x(i-1));
grad_u5(j,i) = (u5(j,i+1)-u5(j,i-1))/(x(i+1)-x(i-1));
end
end
for i=(xstep/5)+1:xstep
if i == xstep/5+1
grad_u1(j,i) = (u1(j,i+1)-u1(j,i))/(x(i+1)-x(i));
grad_u2(j,i) = (u2(j,i+1)-u2(j,i))/(x(i+1)-x(i));
grad_u3(j,i) = (u3(j,i+1)-u3(j,i))/(x(i+1)-x(i));
grad_u4(j,i) = (u4(j,i+1)-u4(j,i))/(x(i+1)-x(i));
grad_u5(j,i) = (u5(j,i+1)-u5(j,i))/(x(i+1)-x(i));
elseif i == xstep
grad_u1(j,i) = (u1(j,i)-u1(j,i-1))/(x(i)-x(i-1));
grad_u2(j,i) = (u2(j,i)-u2(j,i-1))/(x(i)-x(i-1));
grad_u3(j,i) = (u3(j,i)-u3(j,i-1))/(x(i)-x(i-1));
grad_u4(j,i) = (u4(j,i)-u4(j,i-1))/(x(i)-x(i-1));
grad_u5(j,i) = (u5(j,i)-u5(j,i-1))/(x(i)-x(i-1));
else
grad_u1(j,i) = (u1(j,i+1)-u1(j,i-1))/(x(i+1)-x(i-1));
grad_u2(j,i) = (u2(j,i+1)-u2(j,i-1))/(x(i+1)-x(i-1));
grad_u3(j,i) = (u3(j,i+1)-u3(j,i-1))/(x(i+1)-x(i-1));
grad_u4(j,i) = (u4(j,i+1)-u4(j,i-1))/(x(i+1)-x(i-1));
grad_u5(j,i) = (u5(j,i+1)-u5(j,i-1))/(x(i+1)-x(i-1));
end
end
end
for i=1:tstep
for j=1:xstep
v_e = -grad_u1(i,j)*grad_u2(i,j)/u2(i,j)^2; % -g_p*g_n/n^2
flowshear_p(i,j) = 1+ alpha_chi*v_e^2;
flowshear_n(i,j) = 1+ alpha_D*v_e^2;
if abs(grad_u1(i,j)) < abs(data.constant.critgradpressure) %&& abs(grad_u2(i,j)) < abs(data.constant.critgraddensity)
Q(i,j) = -grad_u1(i,j)*chi0;
Q0(i,j) = Q(i,j);
Q1(i,j) = 0;
neo_p(i,j) = chi0*(1+grad_u1(i,j))/(1+grad_u1(i,j));
ano_p(i,j) = 0;
Gam(i,j) = -grad_u2(i,j)*D0;
Gam0(i,j) = Gam(i,j);
Gam1(i,j) = 0;
neo_n(i,j) = D0*(1+grad_u2(i,j))/(1+grad_u2(i,j));
ano_n(i,j) = 0;
elseif abs(grad_u1(i,j)) >= abs(data.constant.critgradpressure) %&& abs(grad_u2(i,j)) < abs(data.constant.critgraddensity)
% Q(i,j) = -grad_u1(i,j)*(chi0+chi1*(abs(grad_u1(i,j))+data.constant.critgradpressure)/flowshear_p(i,j));
Q(i,j) = (chi0*(-grad_u1(i,j)) + chi_ano*(-grad_u1(i,j)+data.constant.critgradpressure)*u3(i,j)/flowshear_p(i,j))*(-grad_u1(i,j));
Q0(i,j) = -grad_u1(i,j)*chi0;
Q1(i,j) = (chi_ano*(-grad_u1(i,j)+data.constant.critgradpressure)*u3(i,j)/flowshear_p(i,j))*(-grad_u1(i,j));
neo_p(i,j) = chi0*(1+grad_u1(i,j))/(1+grad_u1(i,j));
% ano_p(i,j) = chi1*(abs(grad_u1(i,j))+data.constant.critgradpressure)/flowshear_p(i,j);
ano_p(i,j) = chi_ano*(abs(grad_u1(i,j))+data.constant.critgradpressure)*u3(i,j)/flowshear_p(i,j)*(-grad_u1(i,j));
Gam(i,j) = -grad_u2(i,j)*D0;
Gam0(i,j) = Gam(i,j);
Gam1(i,j) = 0;
neo_n(i,j) = D0*(1+grad_u2(i,j))/(1+grad_u2(i,j));
ano_n(i,j) = 0;
% elseif abs(grad_u1(i,j)) < abs(data.constant.critgradpressure) %&& abs(grad_u2(i,j)) >= abs(data.constant.critgraddensity)
% Q(i,j) = -grad_u1(i,j)*chi0;
% Q0(i,j) = Q(i,j);
% Q1(i,j) = 0;
% neo_p(i,j) = chi0*(1+grad_u1(i,j))/(1+grad_u1(i,j));
% ano_p(i,j) = 0;
%Gam(i,j) = -grad_u2(i,j)*(D0+D1*(abs(grad_u2(i,j))+data.constant.critgraddensity)/flowshear_n(i,j));
% Gam(i,j) = -grad_u2(i,j)*(D0 + D_ano*(-grad_u2(i,j)+data.constant.critgraddensity)*u3(i,j)/flowshear_n(i,j));
% Gam0(i,j) = -grad_u2(i,j)*D0;
% Gam1(i,j) = -grad_u2(i,j)*(D_ano*(-grad_u2(i,j)+data.constant.critgraddensity)*u3(i,j)/flowshear_n(i,j));
% neo_n(i,j) = D0*(1+grad_u2(i,j))/(1+grad_u2(i,j));
% ano_n(i,j) = D_ano*(abs(grad_u2(i,j))+data.constant.critgraddensity)/flowshear_n(i,j);
else
% Q(i,j) = -grad_u1(i,j)*(chi0+chi1*(abs(grad_u1(i,j))+data.constant.critgradpressure)/flowshear_p(i,j));
Q(i,j) = (chi0*(-grad_u1(i,j)) + chi_ano*(-grad_u1(i,j)+data.constant.critgradpressure)*u3(i,j)/flowshear_p(i,j))*(-grad_u1(i,j));
Q0(i,j) = -grad_u1(i,j)*chi0;
%Q1(i,j) = -grad_u1(i,j)*chi1*(abs(grad_u1(i,j))+data.constant.critgradpressure)/flowshear_p(i,j);
Q1(i,j) = (chi_ano*(-grad_u1(i,j)+data.constant.critgradpressure)*u3(i,j)/flowshear_p(i,j))*(-grad_u1(i,j));
neo_p(i,j) = chi0*(1+grad_u1(i,j))/(1+grad_u1(i,j));
% ano_p(i,j) = chi1*(abs(grad_u1(i,j))+data.constant.critgradpressure)/flowshear_p(i,j);
ano_p(i,j) = chi_ano*(abs(grad_u1(i,j))+data.constant.critgradpressure)*u3(i,j)/flowshear_p(i,j)*(-grad_u1(i,j));
%Gam(i,j) = -grad_u2(i,j)*(D0+D1*(abs(grad_u2(i,j))+data.constant.critgraddensity)/flowshear_n(i,j));
Gam(i,j) = -grad_u2(i,j)*(D0 + D_ano*(-grad_u2(i,j)+data.constant.critgraddensity)*u3(i,j)/flowshear_n(i,j));
Gam0(i,j) = -grad_u2(i,j)*D0;
Gam1(i,j) = -grad_u2(i,j)*D_ano*(-grad_u2(i,j)+data.constant.critgraddensity)*u3(i,j)/flowshear_n(i,j);
neo_n(i,j) = D0*(1+grad_u2(i,j))/(1+grad_u2(i,j));
ano_n(i,j) = D_ano*(abs(grad_u2(i,j))+data.constant.critgraddensity)/flowshear_n(i,j);
end
end
end
%for i=1:xstep
% a5(i) = sol_dirac(i) ;
% data.variable.a5_term(i) = a5(i);
%end
%curve_u = gradient(grad_u,(x(2)-x(1)));
for j=1:tstep
for i = 1:xstep/5
if i == 1
curve_u1(j,i) = (grad_u1(j,2)-grad_u1(j,1))/(x(2)-x(1));
curve_u2(j,i) = (grad_u2(j,2)-grad_u2(j,1))/(x(2)-x(1));
curve_u3(j,i) = (grad_u3(j,2)-grad_u3(j,1))/(x(2)-x(1));
curve_u4(j,i) = (grad_u4(j,2)-grad_u4(j,1))/(x(2)-x(1));
curve_u5(j,i) = (grad_u5(j,2)-grad_u5(j,1))/(x(2)-x(1));
elseif i == xstep/5
curve_u1(j,i) = (grad_u1(j,i)-grad_u1(j,i-1))/(x(i)-x(i-1));
curve_u2(j,i) = (grad_u2(j,i)-grad_u2(j,i-1))/(x(i)-x(i-1));
curve_u3(j,i) = (grad_u3(j,i)-grad_u3(j,i-1))/(x(i)-x(i-1));
curve_u4(j,i) = (grad_u4(j,i)-grad_u4(j,i-1))/(x(i)-x(i-1));
curve_u5(j,i) = (grad_u5(j,i)-grad_u5(j,i-1))/(x(i)-x(i-1));
else
curve_u1(j,i) = (grad_u1(j,i+1)-grad_u1(j,i-1))/(x(i+1)-x(i-1));
curve_u2(j,i) = (grad_u2(j,i+1)-grad_u2(j,i-1))/(x(i+1)-x(i-1));
curve_u3(j,i) = (grad_u3(j,i+1)-grad_u3(j,i-1))/(x(i+1)-x(i-1));
curve_u4(j,i) = (grad_u4(j,i+1)-grad_u4(j,i-1))/(x(i+1)-x(i-1));
curve_u5(j,i) = (grad_u5(j,i+1)-grad_u5(j,i-1))/(x(i+1)-x(i-1));
end
end
for i=(xstep/5)+1:xstep
if i == xstep/5+1
curve_u1(j,i) = (grad_u1(j,i+1)-grad_u1(j,i))/(x(i+1)-x(i));
curve_u2(j,i) = (grad_u2(j,i+1)-grad_u2(j,i))/(x(i+1)-x(i));
curve_u3(j,i) = (grad_u3(j,i+1)-grad_u3(j,i))/(x(i+1)-x(i));
curve_u4(j,i) = (grad_u4(j,i+1)-grad_u4(j,i))/(x(i+1)-x(i));
curve_u5(j,i) = (grad_u5(j,i+1)-grad_u5(j,i))/(x(i+1)-x(i));
elseif i == xstep
curve_u1(j,i) = (grad_u1(j,i)-grad_u1(j,i-1))/(x(i)-x(i-1));
curve_u2(j,i) = (grad_u2(j,i)-grad_u2(j,i-1))/(x(i)-x(i-1));
curve_u3(j,i) = (grad_u3(j,i)-grad_u3(j,i-1))/(x(i)-x(i-1));
curve_u4(j,i) = (grad_u4(j,i)-grad_u4(j,i-1))/(x(i)-x(i-1));
curve_u5(j,i) = (grad_u5(j,i)-grad_u5(j,i-1))/(x(i)-x(i-1));
else
curve_u1(j,i) = (grad_u1(j,i+1)-grad_u1(j,i-1))/(x(i+1)-x(i-1));
curve_u2(j,i) = (grad_u2(j,i+1)-grad_u2(j,i-1))/(x(i+1)-x(i-1));
curve_u3(j,i) = (grad_u3(j,i+1)-grad_u3(j,i-1))/(x(i+1)-x(i-1));
curve_u4(j,i) = (grad_u4(j,i+1)-grad_u4(j,i-1))/(x(i+1)-x(i-1));
curve_u5(j,i) = (grad_u5(j,i+1)-grad_u5(j,i-1))/(x(i+1)-x(i-1));
end
end
end
%to save parameters and variables
data.constant.chi0 = chi0;
data.constant.D0 = D0;
data.constant.chi1 = chi1;
data.constant.D1 = D1;
data.constant.alphachi = alpha_chi;
data.constant.alphaD = alpha_D;
data.constant.H0 = H0;
data.constant.S0 = S0;
data.variable.pressure = u1;
data.variable.density = u2;
data.variable.turbulence = u3;
data.variable.zonalflow = u4;
data.variable.meanflow = u5;
data.variable.gradpressure = grad_u1;
data.variable.graddensity = grad_u2;
data.variable.gradintensity = grad_u3;
data.variable.curvepressure = curve_u1;
data.variable.curvedensity = curve_u2;
data.variable.curveturbulence = curve_u3;
data.variable.x = x;
data.variable.t = t;
data.variable.Q = Q;
data.variable.Gamma = Gam;
data.variable.Q0 = Q0;
data.variable.Gamma0 = Gam0;
data.variable.neo_P = neo_p;
data.variable.neo_n = neo_n;
data.variable.Q1 = Q1;
data.variable.Gam1 = Gam1;
data.variable.ano_p = ano_p;
data.variable.ano_n = ano_n;
data.variable.heatsource = H;
data.variable.particlesource = S;
data.variable.wexb_p = flowshear_p;
data.variable.wexb_n = flowshear_n;
%data.predicted.Q = Q_pre;
%data.predicted.Gam = Gam_pre;
%data.predicted.gradientp = pp_pre;
%data.predicted.gradientn = nn_pre;
data.control.xgrid = xstep;
data.control.tgrid = tstep;
data.control.xmin = xmin;
data.control.xmax = xmax;
data.control.tmin = tmin;
data.control.tmax = tmax;
figure;
subplot(2,1,1,'FontSize',22)
plot(x,u1(end,:),'.');
axis ([0 1 0 20]);
ylabel('p')
xlabel('r/a')
%set(gca,'XTick',[0,0.2,0.4,0.6,0.8])
head = strcat('Time = ',num2str(t(end)),' s');
title(head);
grid on
subplot(2,1,2,'FontSize',22)
plot(x,grad_u1(end,:),'.');
axis([0 1 -60 0]);
ylabel('p\prime')
xlabel('r/a')
grid on
figure;
subplot(1,2,1,'FontSize',22)
plot(x,curve_u1(end,:),'.');
axis ([0 1 -1400 0]);
ylabel('p\prime\prime')
xlabel('r/a')
%set(gca,'XTick',[0,0.2,0.4,0.6,0.8])
head = strcat('Time = ',num2str(t(end)),' s');
title(head);
grid on
subplot(1,2,2,'FontSize',22)
plot(abs(grad_u1(end,:)),Q(end,:),'.');
%plot(pp_pre,Q_pre,'color','blue','LineWidth',3);
%hold on
%axis ([0 60 0 60]);
ylabel('Q')
xlabel('|p\prime|')
%title('Bifurcation diagram')
grid on
figure;
subplot(2,1,1,'FontSize',22)
plot(x,u2(end,:),'.');
axis ([0 1 0 20]);
ylabel('n')
xlabel('r/a')
%set(gca,'XTick',[0,0.2,0.4,0.6,0.8])
head = strcat('Time = ',num2str(t(end)),' s');
title(head);
grid on
subplot(2,1,2,'FontSize',22)
plot(x,grad_u2(end,:),'.');
axis([0 1 -60 0]);
ylabel('n\prime')
xlabel('r/a')
grid on
figure;
subplot(1,2,1,'FontSize',22)
plot(x,curve_u2(end,:),'.');
axis ([0 1 -1400 0]);
ylabel('n\prime\prime')
xlabel('r/a')
%set(gca,'XTick',[0,0.2,0.4,0.6,0.8])
head = strcat('Time = ',num2str(t(end)),' s');
title(head);
grid on
subplot(1,2,2,'FontSize',22)
plot(abs(grad_u2(end,:)),Gam(end,:),'.');
%plot(nn_pre,Gam_pre,'color','blue','LineWidth',3);
%hold on
%axis ([0 60 0 60]);
ylabel('\Gamma')
xlabel('|n\prime|')
%title('Bifurcation diagram')
grid on
figure;
subplot(2,1,1,'FontSize',22)
plot(x,u3(end,:),'.');
axis ([0 1 0 20]);
ylabel('I')
xlabel('r/a')
%set(gca,'XTick',[0,0.2,0.4,0.6,0.8])
head = strcat('Time = ',num2str(t(end)),' s');
title(head);
grid on
subplot(2,1,2,'FontSize',22)
plot(x,grad_u3(end,:),'.');
axis([0 1 -60 0]);
ylabel('I\prime')
xlabel('r/a')
grid on
figure;
subplot(1,2,1,'FontSize',22)
plot(x,curve_u3(end,:),'.');
axis ([0 1 -1400 0]);
ylabel('I\prime\prime')
xlabel('r/a')
%set(gca,'XTick',[0,0.2,0.4,0.6,0.8])
head = strcat('Time = ',num2str(t(end)),' s');
title(head);
grid on
subplot(1,2,2,'FontSize',22)
plot(abs(grad_u3(end,:)),Q(end,:),'.');
%plot(pp_pre,Q_pre,'color','blue','LineWidth',3);
%hold on
%axis ([0 60 0 60]);
ylabel('Q')
xlabel('|I\prime|')
%title('Bifurcation diagram')
grid on
end
% --------------------------------------------------------------
function [c,f,s] = pdex2pde(x,t,u,DuDx)
global chi0;
global D0;
global chi1;
global D1;
global H0;
global S0;
global data;
global track;
global xstep;
global tstep;
global alpha_chi;
global alpha_D;
global chi_growth; % total growth rate
global length;
global theta_heaviside1;
global v_e;
global chi_ano;
global D_ano;
global sigma_turb;
global gamma_nonlin;
global alpha_zero;
global alpha_v;
global rho_0;
global gamma_damp;
global alpha_nonlin;
global dCoeff;
%lf FFf F Fb vfDdength = 0.01 or 1
length=1;
c = [1;1;1;1;1];
v_e = -DuDx(1)*DuDx(2)/u(2)^2; % -g_p*g_n/n^2
flowshear_p = 1+ alpha_chi*v_e^2;
flowshear_n = 1+ alpha_D*v_e^2;
term1 = abs(DuDx(1))-data.constant.critgradpressure;
% Implementing Heaviside function for H-mode in p, n and I equations
if term1 > 0
theta_heaviside1=1;
else
theta_heaviside1=0;
end
%Turbulence intensity Equation for random walk
s_1 = alpha_zero*u(4)*u(3)/(1+rho_0*u(5))-gamma_damp*u(4);
s_2 = mean(v_e)^2;
s_3 = chi_growth*(term1*theta_heaviside1-gamma_nonlin*u(3)^alpha_nonlin-alpha_zero*u(4)-alpha_v*u(5))*u(3) ;
s = [(H0)*exp(-100*x^2/length)+H0/2; (S0)*exp(-100*(x-0.9)^2/length)+S0/2; s_1;s_2;s_3];
f = [chi0+chi1*u(3)/flowshear_p ; D0+D1*u(3)/flowshear_n ; 0; 0; dCoeff].*DuDx; % flux term for random walk model
%data.function.c(track) = c;
%data.function.f(track) = f;
if track <= xstep
data.function.heatsource(track) = s(1);
data.function.particlesource(track) = s(2);
data.function.turbulencesource(track)= s(3);
data.function.x(track) = x;
track = track+1;
end
end
% --------------------------------------------------------------
function u0 = pdex2ic(x)
%u0 = [eps; eps; eps];
%u0 = [0.01; 0.01; 0.1*exp(-100*(x-1)^2)];
u0 = [0.4*(1-x^2); 0.4*(1-x^2); 0.5*exp(-100*(x-1)^2);0.4*(1-x^2);0.4*(1-x^2)];
end
% --------------------------------------------------------------
function [pl,ql,pr,qr] = pdex2bc(xl,ul,xr,ur,t)
pl = [0; 0; 0; 0; 0];
ql = [1; 1; 1; 1; 1];
pr = [ur(1); ur(2); ur(3); ur(4);0];
qr = [0.01; 0.1; 1;1;1];
% qr = [0; 0.1; 1];
%qr = [0; eps; 1];
%---------------------------------------------------------------
end
Torsten
Torsten 2023-10-18
编辑：Torsten 2023-10-18
Your code doesn't produce syntax errors up to the point where you call "pdepe". The solver starts integrating, but stops after 1.855e-5 s because it cannot continue. This is usually an indication that either the equations are hard to solve or you made some errors in the setup.
I suggest you include the mathematical formulation of the equations you are trying to solve with "pdepe".

请先登录，再进行评论。

类别

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

标签

产品


版本

R2023b

Community Treasure Hunt

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

Start Hunting!

Translated by