How to fix the error in pdepe?

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()

Warning: The value of local variables may have been changed to match the globals. Future versions of MATLAB will require that you declare a variable to be global before you use that variable.

Warning: The value of local variables may have been changed to match the globals. Future versions of MATLAB will require that you declare a variable to be global before you use that variable.

     0

   2.2204e-16

ans = 1×2

3 1

ans = 1×2

3 1

ans = 1×2

1002 100

Error using pdepe
Unexpected output of PDEFUN. For this problem PDEFUN must return three column vectors of length 3.

Error in solution>pde2fshear_v4Perturbed_nonlinear (line 152)
sol = pdepe(m,@pdex2pde,@pdex2ic,@pdex2bc,x,t);

% 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

Answer 1

[Torsten]

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.

[Rahul]

Hi Torsten,

Thanks a lot for your valuable advice.

with regards,

rc