How to fix the error in pdepe?
2 visualizzazioni (ultimi 30 giorni)
Mostra commenti meno recenti
Hi,
In the attached code, it throws an error which says
Error using ^
Incorrect dimensions for raising a matrix to a power. Check that the matrix is square and the power is a scalar. To
operate on each element of the matrix individually, use POWER (.^) for elementwise power.
Now, when I try to do the correction by putting .^ in place of ^ for vector multiplication of vector x, it again throws an error which says
Error using pdepe
Unexpected output of PDEFUN. For this problem PDEFUN must return three column vectors of length 3.
I need your advice.
pde2fshear_v4Perturbed_nonlinear()
% solve 3-F bifurcaiton model
function pde2fshear_v4Perturbed_nonlinear
global x;
global X;
global t;
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 chi_growth;
global lambda_suppress;
global v_e;
global chi_ano;
global D_ano;
global gamma_nonlin;
global alpha_nonlin;
global group_vel;
global drift_vel;
global intensity_diff;
global drift_velFluct;
global D_0;
global safetyFact;
global grad_safetyFact;
global elec_diaVel;
global grad_length;
global mag_shear;
global theta_deg;
global elec_temp;
global T_magFld;
global P_magFld;
global R0;
global I_p;
global jb;
global jd;
global jtot;
global B0;
global xmax;
global xmin;
global STEP_x;
%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 = 27;
S0 = 21;
track = 1;
chi_growth=20;
lambda_suppress=0.05;
chi_ano=10;
D_ano=10;
% For nonlinear model
gamma_nonlin = 50;
alpha_nonlin = 1;
D_0 = 10;
STEP_x=1;
%position and time grids information
xstep = 100;
tstep = 1000;
xmin = 0;
%xmin = eps;
xmax = 1;
tmin = 0;
%tmin = eps;
tmax = 50;
%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);
curve_u1 = zeros(tstep,xstep);
curve_u2 = zeros(tstep,xstep);
curve_u3 = 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);
safetyFact=zeros(tstep,xstep);
elec_diaVel=zeros(tstep,xstep);
grad_length=zeros(tstep,xstep);
mag_shear=zeros(tstep,xstep);
theta_deg=zeros(tstep,xstep);
elec_temp=zeros(tstep,xstep);
T_magFld=zeros(tstep,xstep);
P_magFld=zeros(tstep,xstep);
R0=zeros(tstep,xstep);
I_p=zeros(tstep,xstep);
%jb=zeros(tstep,xstep);
B0=zeros(tstep,xstep);
%define first inner half of the plasma
%x = linspace(xmin,xmax/2,xstep/5);
x = linspace(xmin,xmax,xstep);
X = linspace(xmin,xmax,xstep);
Xi=xmin;
Xf=xmax;
XRange=[Xi,Xf];
idl = X >= XRange(1) & X <= XRange(2);
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
u1 = sol(:,:,1);
u2 = sol(:,:,2);
u3 = sol(:,:,3);
%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));
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));
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));
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));
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));
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));
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;
intensity_diff(i,j) = D_0*u3(i,j).^alpha_nonlin;
drift_velFluct(i,j)=D_0*alpha_nonlin*grad_u3(i,j).^(alpha_nonlin-1)*grad_u3(i,j);
drift_vel(i,j) = group_vel + drift_velFluct(i,j);
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)) + 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) = 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;
intensity_diff(i,j) = D_0*u3(i,j).^alpha_nonlin;
drift_velFluct(i,j)=D_0*alpha_nonlin*grad_u3(i,j).^(alpha_nonlin-1)*grad_u3(i,j);
drift_vel(i,j) = group_vel + drift_velFluct(i,j);
else
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) = 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 + 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);
intensity_diff(i,j) = D_0*u3(i,j).^alpha_nonlin;
drift_velFluct(i,j)=D_0*alpha_nonlin*grad_u3(i,j).^(alpha_nonlin-1)*grad_u3(i,j);
drift_vel(i,j) = group_vel + drift_velFluct(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));
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));
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));
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));
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));
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));
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.intensity_diff=intensity_diff;
data.variable.drift_velFluct=drift_velFluct;
data.variable.drift_vel=drift_vel;
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.control.xgrid = xstep;
data.control.tgrid = tstep;
data.control.xmin = xmin;
data.control.xmax = xmax;
data.control.tmin = tmin;
data.control.tmax = tmax;
end
% --------------------------------------------------------------
function [c,f,s] = pdex2pde(x,t,u,DuDx)
global x;
global t;
global chi0;
global D0;
global chi1;
global D1;
global H0;
global S0;
global data;
global xstep;
global tstep;
global alpha_chi;
global alpha_D;
global chi_growth; % total growth rate
global length;
global theta_heaviside1;
global lambda_suppress;
global v_e;
global gamma_nonlin;
global alpha_nonlin;
global group_vel;
global drift_vel;
global intensity_diff;
global drift_velFluct;
global D_0;
global safetyFact;
global grad_safetyFact;
global elec_diaVel;
global grad_length;
global mag_shear;
global theta_deg;
global elec_temp;
global T_magFld;
global P_magFld;
global R0;
global I_p;
global jb;
global jd;
global jtot;
global B0;
global X;
global xmax;
global xmin;
global STEP_x;
%lf FFf F Fb vfDdength = 0.01 or 1
length=1;
jb0=1;
jd0=1;
x0=1;
grad_length=DuDx(2);
R0=linspace(1,5,100);
B0=1;
c = [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;
drift_velFluct = D_0*DuDx(3)^alpha_nonlin;
intensity_diff = D_0*u(3)^alpha_nonlin;
% Implementing Heaviside function for H-mode in p, n and I equations
if term1 > 0
theta_heaviside1=1;
else
theta_heaviside1=0;
end
% Determining group velocity of turbulence intensity
if STEP_x > xstep-1
STEP_x=1;
end
if STEP_x <= xstep-1
X(STEP_x)=x(STEP_x);
end
disp(X(STEP_x))
num=15;
if STEP_x==1
jb=jb0*DuDx(1);
jd=jd0.*(1-(X(STEP_x)-x0).^2).^num;
jtot=jd+jb(STEP_x);
I_p=eps;
elseif STEP_x==2
jd=j0.*(1-(X(STEP_x)-x0).^2).^num;
jb(STEP_x)=jb0*DuDx(1);
jtot=jd+jb(STEP_x);
I_p=trapz(X(1:STEP_x),X(1:STEP_x).*jb(1:STEP_x) + X(1:STEP_x).*j0.*(1-(X(1:STEP_x)-x0).^2).^num);
elseif STEP_x==3
jd=j0.*(1-(X(STEP_x)-x0).^2).^num;
jb(STEP_x)=jb0*DuDx(1);
jtot=jd+jb(STEP_x);
I_p=trapz(X(1:STEP_x),X(1:STEP_x).*(jb(1:STEP_x)+j0.*(1-(X(1:STEP_x)-x0).^2).^num));
else
jd=j0.*(1-(X(STEP_x)-x0).^2).^nu;
jb(STEP_x)=jb0*DuDx(1);
jtot=jd+jb(STEP_x);
I_p=trapz(X(1:STEP_x),X(1:STEP_x).*(jb(1:STEP_x)+j0.*(1-(X(1:STEP_x)-x0).^2).^num));
end
T_magFld=B0./R0;
P_magFld=B0.*I_p./x;
safetyFact=(x./R0).*T_magFld./P_magFld;
if STEP_x == 1
grad_safetyFact(STEP_x) = (safetyFact(2)-safetyFact(1))/(x(2)-x(1));
elseif STEP_x == xstep/5
grad_safetyFact(STEP_x) = (safetyFact(STEP_x)-safetyFact(STEP_x-1))/(x(STEP_x)-x(STEP_x-1));
else
grad_safetyFact(STEP_x) = (safetyFact(STEP_x+1)-safetyFact(STEP_x-1))/(x(STEP_x+1)-x(STEP_x-1));
end
mag_shear=(x/safetyFact)*grad_safetyFact;
%disp (mag_shear)
elec_diaVel = -(elec_temp/(1.6*10^(-19).*T_magFld)).*grad_length;
%disp(elec_diaVel)
group_vel = - (elec_diaVel.*(2*grad_length/mag_shear.*R0).*sin(theta_deg));
%disp(group_vel)
disp(I_p)
%Turbulence intensity Equation for nonlinear turbulence intensity
s3 = (chi_growth*(term1*theta_heaviside1-lambda_suppress*v_e^2)-gamma_nonlin*u(3)^alpha_nonlin)*u(3)+group_vel*u(3)-(D_0*(1+alpha_nonlin)*alpha_nonlin*u(3)^(alpha_nonlin-1)*u(3)^2);
s = [(H0)*exp(-100*x.^2/length)+H0/2; (S0)*exp(-100*(x-0.9).^2/length)+S0/2;s3];
f = [chi0+chi1*u(3)/flowshear_p ; D0+D1*u(3)/flowshear_n ;-(intensity_diff-D_0*(1+alpha_nonlin)*u(3)^alpha_nonlin)].*DuDx; % flux term for nonlinear model
size(c)
size(f)
size(s)
end
% --------------------------------------------------------------
function u0 = pdex2ic(x)
%u0 = [eps; eps; eps];
%u0 = [0.01; 0.01;0.1*exp(-100*(x-1)^2)];
u0 = [0.1*(1-x^2); 0.1*(1-x^2); 0.5*exp(-100*(x-1)^2)];
end
% --------------------------------------------------------------
function [pl,ql,pr,qr] = pdex2bc(xl,ul,xr,ur,t)
pl = [0; 0; 0];
ql = [1; 1; 1];
pr = [ur(1); ur(2);0];
qr = [0.01; 0.1; 1];
%---------------------------------------------------------------
end
with regards,
rc
0 Commenti
Risposta accettata
Torsten
il 4 Feb 2024
Modificato: Torsten
il 4 Feb 2024
I included the commands
size(c)
size(f)
size(s)
at the end of your function pdex2pde.
You can see that MATLAB doesn't lie: c and f are of size 3x1 as required, but s is of size 1002x100.
Since nobody will be able to understand what you are doing in the code, this is just to show your error, not to supply an answer.
Più risposte (0)
Vedere anche
Categorie
Scopri di più su Numerical Integration and Differential Equations in Help Center e File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!