I get the error but I dont know how to fix it even try:(

Question

sevgul demir 2021-11-25

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1595599-i-get-the-error-but-i-dont-know-how-to-fix-it-even-try

评论： Dyuman Joshi 2021-11-26

error.mlx

% Commented Matlabfi code for the resolution of the sixth case stdUY
clear all
close all
% Computation precision
tol=1e-5;
% x, y and z dimension of the 3D domain. for this case, its a cubic
% domain
Ax=1; Ay=1; Az=1;
% Temperature of the hot (Tamx) & cold (Tmin) sidewalls in Kelvin
Tmin=400; Tmax=1200;
% Inputs for the thermal conductivity as a function of temperature
lx=420.75; ly=420.75; lz=420.75; ax=-(1.627878)*(1e-4);
ay=-(1.627878)*(1e-4); az=-(1.627878)*(1e-4); H1=15;
H3=15;
% number of collocation points in the x, y and z directions
Nx=12;Ny=12;Nz=12;
% 1D First and second order derivatives and identity matrixes (see
% N. L. Trefethen, Spectral Methods in MATLAB, Philadelphia:
% SIAM, 2000)
x = cos (pi*(0:Nx)/Nx);
c = [2; ones(Nx-1,1); 2].*(-1).^(0:Nx);
X = repmat(x,1,Nx+1);
dX = X-X;
Dx = (c*(1./c))./(Dx+(eye(Nx+1)));
Dx = Dx - diag(sum(Dx));
D1_x = Dx; D2_x = Dx^2; I_x = speye(Nx+1);
y = cos (pi*(0:Ny)/Ny);
c = [2; ones(Ny-1,1); 2].*(-1).^(0:Ny);
Y = repmat(y,1,Ny+1);
dY = Y-Y;
Dy = (c*(1./c))./(dY+(eye(Nx+1)));
Dy = Dy - diag(sum(Dy));
D1_y = Dy; D2_y = Dy^2; I_y = speye(Ny+1);
z = cos (pi*(0:Nz)/Nz);
c = [2; ones(Nz-1,1); 2].*(-1).^(0:Nz);
Z = repmat(z,1,Nz+1);
dZ = Z-Z;
Dz = (c*(1./c))./(dZ+(eye(Nz+1)));
Dz = Dz - diag(sum(Dz));
D1_z = Dz; D2_z = Dz^2; I_z = speye(Nz+1);
% Adjust the 1D mesh collocation points
x=Ax*x; y=Ay*y; z=Az*z;
% 3D meshgrid
[xx,yy,zz] = meshgrid(x,y,z);
% Vectorize the mesh points
xxv = xx(:); yyv = yy(:); zzv = zz(:);
% 3D partial derivatives with respect to x, y and z directions
% 3D-First order partial derivative with respect to x: D/Dx
D1x=(1/(Ax))*kron(kron(I_z,D1_x),I_y);
% 3D-First order partial derivative with respect to y: D/Dy
D1y = (1/(Ay))*kron(kron(I_x,I_z),D1_y);
% 3D-First order partial derivative with respect to z: D/Dz
D1z= (1/(Az))*kron(kron(D1_z,I_x),I_y);
% 3D-Second order partial derivative with respect to x: D2/Dx2
D2x= (1/(Ax^2))*kron(kron(I_z,D2_x),I_y);
% 3D-Second order partial derivative with respect to y: D2/Dy2
D2y= (1/(Ay^2))*kron(kron(I_x,I_z),D2_y);
% 3D-Second order partial derivative with respect to z: D2/Dz2
D2z= (1/(Az^2))*kron(kron(D2_z,I_x),I_y);
% Construction of the 3D Laplacian
L3=D2x+D2y+D2z;
% Identity and zeros matrixes
Ixyz=speye(N1x*N1y*N1z);
Oxyz=0*Ixyz;
% Initialization of the temperature (between Tmin=400 &Tmax=1200 K)
k1 = (Tmax - Tmin)/2;
k2 = (Tmax + Tmin)/2;
T = k2 + k1*(yyv/Ay);
% Location of the boundaries
b = find(abs(xx)==Ax |abs(yy)==Ay | abs(zz)==Az);
bx = find(abs(xx)==Ax);
bx_1 = find(xx==-Ax);
bx1 = find(xx==Ax);
by = find(abs(yy)==Ay);
by_1 = find(yy==-Ay & abs(xx)<Ax & abs(zz)<Az);
by1 = find(yy==Ay & abs(xx)<Ax & abs(zz)<Az);
bz = find(abs(zz)==Az);
bz_1 = find(zz==-Az);
bz1 = find(zz==Az);
% Start of the fixed point method iterations
% initialize the difference between the temperature values between two
% successive iterations
dif_it = 10; it = 0;
while dif_it > tol
    % 3D operator for the thermal conduction problem with variable thermal
    % conductivity: Ltemp
    lamx=lx*(ones(N1x*N1y*N1z,1)+(ax*T));
    lamy=ly*(ones(N1x*N1y*N1z,1)+(ay*T));
    lamz=lz*(ones(N1x*N1y*N1z,1)+(az*T));
    Ltemp=lx*ax*diag(D1x*T)*D1x+diag(lamx)*D2x+ly*ay*diag(D1y*T)*D1y+diag(lamy)*D2y+lz*az*diag(D1z*T)*D1z+diag(lamz)*D2z;
    D1xCL=diag(lamx)*D1x;
    D1yCL=diag(lamy)*D1y;
    D1zCL=diag(lamz)*D1z;
    % Enforce boudary conditions on the Ltemp operator
    Ltemp(bx1,:) = Oxyz(bx1,:); Ltemp(bx1,:) = D1xCL(bx1,:)+H1*Ixyz(bx1,:);
    Ltemp(bx_1,:) = Oxyz(bx_1,:); Ltemp(bx_1,:) = D1xCL(bx_1,:)-H1*Ixyz(bx_1,:);
    Ltemp(by_1,:) = Ixyz(by_1,:);
    Ltemp(by1,:) = Ixyz(by1,:);
    Ltemp(bz1,:) = Oxyz(bz1,:); Ltemp(bz1,:) = D1zCL(bz1,:)+H3*Ixyz(bz1,:);
    Ltemp(bz_1,:) = Oxyz(bz_1,:); Ltemp(bz_1,:) = D1zCL(bz_1,:)-H3*Ixyz(bz_1,:);
    % 3D operator for the thermal conduction problem with constant thermal
    % conduction: Lcond
    Lcond=L3;
    % Enforce boudary conditions on the Lcond operator
    Lcond(bx1,:) = Oxyz(bx1,:); Lcond(bx1,:) =lx*D1x(bx1,:)+H1*Ixyz(bx1,:);
    Lcond(bx_1,:) = Oxyz(bx_1,:);
    Lcond(bx_1,:) = lx*D1x(bx_1,:)-H1*Ixyz(bx_1,:);
    Lcond(by_1,:) = Ixyz(by_1,:);
    Lcond(by1,:) = Ixyz(by1,:);
    Lcond(bz1,:) = Oxyz(bz1,:); Lcond(bz1,:) =lx*D1z(bz1,:)+H3*Ixyz(bz1,:);
    Lcond(bz_1,:) = Oxyz(bz_1,:);
    Lcond(bz_1,:) = lx*D1z(bz_1,:)-H3*Ixyz(bz_1,:);
    % Right hand term
    rhstemp=zeros(N1x*N1y*N1z,1);
    % Enforce boudary conditions on the righthand term
    rhstemp(bx1) = 300*H1;
    rhstemp(bx_1) = -300*H1;
    rhstemp(by_1) =400;
    rhstemp(by1) = 1200;
    rhstemp(bz1) = 300*H3;
    rhstemp(bz_1) = -300*H3;
    % Solution with variable thermal conductivity
    Tnew=Ltemp\rhstemp;
    % Solution with constant thermal conductivity
    Tc=Lcond\rhstemp;
    % Convergence criterion
    dif_T=Tnew-T;
    dif_it = abs(norm(dif_T,inf));
    % Update the temperature value
    T=Tnew;
    % Update the itration value
    it = it+1;
end