I wrote the below code,but I’m getting the following erroes, how…
QuestionI wrote the below code,but I’m getting the following erroes, how…I wrote the below code,but I’m getting the following erroes, how can that be fixed? Error in HW7>createMesh3D (line 394)G=reshape(1:(Nx+2)*(Ny+2)*(Nz+2), Nx+2, Ny+2, Nz+2); Error in HW7>createMesh3D (line 408)MS=createMesh3D(3,[Nx,Ny, Nz],CellSize,CellLocation,FaceLocation, c(:),[e1(:); e2(:); e3(:)]); Error in HW7 (line 350)m = createMesh3D(N-2,N-2,N-2,Lx,Ly,Lz); % create the mesh HW7394 G=reshape(1:(Nx+2)*(Ny+2)*(Nz+2), Nx+2, Ny+2, Nz+2); Code: .rtcContent { padding: 30px; } .lineNode {font-size: 10pt; font-family: Menlo, Monaco, Consolas, “Courier New”, monospace; font-style: normal; font-weight: normal; }%CP7%Written by Viridiana Salazarclear variablesclose all%% 1. Space discretizationLx = 1.0; dx =0.125; N=Lx/dx+1; x=0:dx:Lx;Ly = 1.0; dy =0.025; y=0:dy:Ly;Lz = 1.0; dz =0.025; z=0:dz:Lz;%% 2. Time discretizationtf =0.5; dt =0.001; M=tf/dt+1; t= 0:dt:tf;%% ConstantsMu=0.5; r=Mu*dt/(dx)^2; q=-0.5*dt/dx;%% Analytical Solution (We’re gonna take from here the initial and% boundary conditions)UA = zeros(N,N,N,M);for n=1:M %timefor i=1:N %xfor j=1:N %yfor k=1:N %zUA(i,j,k,n)=(1+exp(x(i)/(2*Mu)+y(j)/(2*Mu)+z(k)/(2*Mu)-3*t(n)/(4*Mu)))^(-1);endendendend%% Boundary ConditionsU = zeros(N,N,N,M);% For x=1 and x=NU(1,:,:,:)=UA(1,:,:,:); U(N,:,:,:)=UA(N,:,:,:);%For y=1 and y=NU(:,1,:,:)=UA(:,1,:,:); U(:,N,:,:)=UA(:,N,:,:);%For z=1 and z=NU(:,:,1,:)=UA(:,:,1,:); U(:,:,N,:)=UA(:,:,N,:);%% Initial conditions% For t=0U(:,:,:,1)=UA(:,:,:,1);%%for n=1:M-1 %Time loop% Initial guest at level n+1for i=1:N %xfor j=1:N %yfor k=1:N %zif i~=1 && j~=1 && k~=1 && i~=N && j~=N && k~=N %At the boundaries 1 and N, U is knownU(i,j,k,n+1)=U(i,j,k,n);endendendenditeration=0;EPS=1;while EPS>0.0000000001 && iteration<20%Jacobian MatrixJ=zeros(N^3,N^3);% Vector fF=zeros(N^3,1);i=0;for I=1:N^3K=I;% Fix the index i,j,kif mod(I,N*N)==1i=i+1;j=0;endif mod(I,N)==1k=1;j=j+1;elsek=k+1;endif i~=1 && j~=1 && k~=1 && i~=N && j~=N && k~=NJ(K,I-1)=q+r; %DF/U(i,j,k-1)J(K,I)=-3*q-(6*r+1);%DF/U(i,j,k)J(K,I+1)=r; %DF/U(i,j,k+1)J(K,I+N)=r; %DF/U(i,j+1,k)J(K,I-N)=q+r; %DF/U(i,j-1,k)J(K,I+N*N)=r; %DF/U(i+1,j,k)J(K,I-N*N)=q+r; %DF/U(i-1,j,k)F(K,1)=(q*(-3*U(i,j,k,n+1)+U(i-1,j,k,n+1)+U(i,j,k-1,n+1)+U(i,j,k-1,n+1))+(-(6*r+1)*U(i,j,k,n+1)+r*(U(i+1,j,k,n+1)+U(i-1,j,k,n+1)+U(i,j+1,k,n+1)+U(i,j-1,k,n+1)+U(i,j,k+1,n+1)+U(i,j,k-1,n+1)))+U(i,j,k,n));end%% derivatives for i=N and/or j=N and/or k=Nif i==N && j==N && k==NJ(K,I-1)=q+r; %DF/U(i,j,k-1)J(K,I)=-3*q-(6*r+1);%DF/U(i,j,k)J(K,I-N)=q+r; %DF/U(i,j-1,k)J(K,I-N*N)=q+r; %DF/U(i-1,j,k)end%i==N && j~=N && k~=Nif i==N && j~=N && k~=N && j~=1 && k~=1J(K,I-1)=q+r; %DF/U(i,j,k-1)J(K,I)=-3*q-(6*r+1);%DF/U(i,j,k)J(K,I+1)=r; %DF/U(i,j,k+1)J(K,I+N)=r; %DF/U(i,j+1,k)J(K,I-N)=q+r; %DF/U(i,j-1,k)J(K,I-N*N)=q+r; %DF/U(i-1,j,k)endif i==N && j~=N && k~=N && j==1 && k~=1J(K,I-1)=q+r; %DF/U(i,j,k-1)J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)J(K,I+1)=r; %DF/U(i,j,k+1)J(K,I+N)=r; %DF/U(i,j+1,k)J(K,I-N*N)=q+r; %DF/U(i-1,j,k)endif i==N && j~=N && k~=N && j~=1 && k==1J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)J(K,I+1)=r; %DF/U(i,j,k+1)J(K,I+N)=r; %DF/U(i,j+1,k)J(K,I-N)=q+r; %DF/U(i,j-1,k)J(K,I-N*N)=q+r; %DF/U(i-1,j,k)end%i~=N && j==N && k~=Nif i~=N && j==N && k~=N && i~=1 && k~=1J(K,I-1)=q+r; %DF/U(i,j,k-1)J(K,I)=-3*q-(6*r+1);%DF/U(i,j,k)J(K,I+1)=r; %DF/U(i,j,k+1)J(K,I-N)=q+r; %DF/U(i,j-1,k)J(K,I+N*N)=r; %DF/U(i+1,j,k)J(K,I-N*N)=q+r; %DF/U(i-1,j,k)endif i~=N && j==N && k~=N && i==1 && k~=1J(K,I-1)=q*+r; %DF/U(i,j,k-1)J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)J(K,I+1)=r; %DF/U(i,j,k+1)J(K,I-N)=q+r; %DF/U(i,j-1,k)J(K,I+N*N)=r; %DF/U(i+1,j,k)endif i~=N && j==N && k~=N && i~=1 && k==1J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)J(K,I+1)=r; %DF/U(i,j,k+1)J(K,I-N)=q+r; %DF/U(i,j-1,k)J(K,I+N*N)=r; %DF/U(i+1,j,k)J(K,I-N*N)=q+r; %DF/U(i-1,j,k)end%i~=N && j~=N && k==Nif i~=N && j~=N && k==N && i~=1 && j~=1J(K,I-1)=q+r; %DF/U(i,j,k-1)J(K,I)=-3*q-(6*r+1);%DF/U(i,j,k)J(K,I+N)=r; %DF/U(i,j+1,k)J(K,I-N)=q+r; %DF/U(i,j-1,k)J(K,I+N*N)=r; %DF/U(i+1,j,k)J(K,I-N*N)=q+r; %DF/U(i-1,j,k)endif i~=N && j~=N && k==N && i==1 && j~=1J(K,I-1)=q+r; %DF/U(i,j,k-1)J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)J(K,I+N)=r; %DF/U(i,j+1,k)J(K,I-N)=q+r; %DF/U(i,j-1,k)J(K,I+N*N)=r; %DF/U(i+1,j,k)endif i~=N && j~=N && k==N && i~=1 && j==1J(K,I-1)=q+r; %DF/U(i,j,k-1)J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)J(K,I+N)=r; %DF/U(i,j+1,k)J(K,I+N*N)=r; %DF/U(i+1,j,k)J(K,I-N*N)=q+r; %DF/U(i-1,j,k)end% i==N && j==N && k~=Nif i==N && j==N && k~=N && k~=1J(K,I-1)=q+r; %DF/U(i,j,k-1)J(K,I)=-3*q-(6*r+1);%DF/U(i,j,k)J(K,I+1)=r; %DF/U(i,j,k+1)J(K,I-N)=q+r; %DF/U(i,j-1,k)J(K,I-N*N)=q+r; %DF/U(i-1,j,k)endif i==N && j==N && k~=N && k==1J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)J(K,I+1)=r; %DF/U(i,j,k+1)J(K,I-N)=q+r; %DF/U(i,j-1,k)J(K,I-N*N)=q+r; %DF/U(i-1,j,k)end% i==N && j~=N && k==Nif i==N && j~=N && k==N && j~=1J(K,I-1)=q+r; %DF/U(i,j,k-1)J(K,I)=-3*q-(6*r+1);%DF/U(i,j,k)J(K,I+N)=r; %DF/U(i,j+1,k)J(K,I-N)=q+r; %DF/U(i,j-1,k)J(K,I-N*N)=q*+r; %DF/U(i-1,j,k)endif i==N && j~=N && k==N && j==1J(K,I-1)=q+r; %DF/U(i,j,k-1)J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)J(K,I+N)=r; %DF/U(i,j+1,k)J(K,I-N*N)=q+r; %DF/U(i-1,j,k)end% i~=N && j==N && k==Nif i~=N && j==N && k==N && i~=1J(K,I-1)=q+r; %DF/U(i,j,k-1)J(K,I)=-3*q-(6*r+1);%DF/U(i,j,k)J(K,I-N)=q+r; %DF/U(i,j-1,k)J(K,I+N*N)=r; %DF/U(i+1,j,k)J(K,I-N*N)=q+r; %DF/U(i-1,j,k)endif i~=N && j==N && k==N && i==1J(K,I-1)=q+r; %DF/U(i,j,k-1)J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)J(K,I-N)=q+r; %DF/U(i,j-1,k)J(K,I+N*N)=r; %DF/U(i+1,j,k)end %% derivatives for i=1 and/or j=1 and/or k=1if i==1 && j==1 && k==1J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)J(K,I+1)=r; %DF/U(i,j,k+1)J(K,I+N)=r; %DF/U(i,j+1,k)J(K,I+N*N)=r; %DF/U(i+1,j,k)end% i==1 && j~=1 && k~=1 && j~=N && k~=Nif i==1 && j~=1 && k~=1 && j~=N && k~=NJ(K,I-1)=q+r; %DF/U(i,j,k-1)J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)J(K,I+1)=r; %DF/U(i,j,k+1)J(K,I+N)=r; %DF/U(i,j+1,k)J(K,I-N)=q+r; %DF/U(i,j-1,k)J(K,I+N*N)=r; %DF/U(i+1,j,k)endif i==1 && j~=1 && k~=1 && j==N && k~=NJ(K,I-1)=q+r; %DF/U(i,j,k-1)J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)J(K,I+1)=r; %DF/U(i,j,k+1)J(K,I-N)=q+r; %DF/U(i,j-1,k)J(K,I+N*N)=r; %DF/U(i+1,j,k)endif i==1 && j~=1 && k~=1 && j~=N && k==NJ(K,I-1)=q+r; %DF/U(i,j,k-1)J(K,I)=q-(6*r+1); %DF/U(i,j,k)J(K,I+N)=r; %DF/U(i,j+1,k)J(K,I-N)=q+r; %DF/U(i,j-1,k)J(K,I+N*N)=r; %DF/U(i+1,j,k)end% i~=1 && j==1 && k~=1if i~=1 && j==1 && k~=1 && i~=N && k~=NJ(K,I-1)=q+r; %DF/U(i,j,k-1)J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)J(K,I+1)=r; %DF/U(i,j,k+1)J(K,I+N)=r; %DF/U(i,j+1,k)J(K,I+N*N)=r; %DF/U(i+1,j,k)J(K,I-N*N)=q+r; %DF/U(i-1,j,k)endif i~=1 && j==1 && k~=1 && i==N && k~=NJ(K,I-1)=q+r; %DF/U(i,j,k-1)J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)J(K,I+1)=r; %DF/U(i,j,k+1)J(K,I+N)=r; %DF/U(i,j+1,k)J(K,I-N*N)=q+r; %DF/U(i-1,j,k)endif i~=1 && j==1 && k~=1 && i~=N && k==NJ(K,I-1)=q+r; %DF/U(i,j,k-1)J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)J(K,I+N)=r; %DF/U(i,j+1,k)J(K,I+N*N)=r; %DF/U(i+1,j,k)J(K,I-N*N)=q+r; %DF/U(i-1,j,k)end%i~=1 && j~=1 && k==1if i~=1 && j~=1 && k==1 && i~=N && j~=NJ(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)J(K,I+1)=r; %DF/U(i,j,k+1)J(K,I+N)=r; %DF/U(i,j+1,k)J(K,I-N)=q+r; %DF/U(i,j-1,k)J(K,I+N*N)=r; %DF/U(i+1,j,k)J(K,I-N*N)=q+r; %DF/U(i-1,j,k)endif i~=1 && j~=1 && k==1 && i==N && j~=NJ(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)J(K,I+1)=r; %DF/U(i,j,k+1)J(K,I+N)=r; %DF/U(i,j+1,k)J(K,I-N)=q+r; %DF/U(i,j-1,k)J(K,I-N*N)=q+r; %DF/U(i-1,j,k)endif i~=1 && j~=1 && k==1 && i~=N && j==NJ(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)J(K,I+1)=r; %DF/U(i,j,k+1)J(K,I-N)=q+r; %DF/U(i,j-1,k)J(K,I+N*N)=r; %DF/U(i+1,j,k)J(K,I-N*N)=q+r; %DF/U(i-1,j,k) end%i==1 && j==1 && k~=1if i==1 && j==1 && k~=1 && k~=NJ(K,I-1)=q+r; %DF/U(i,j,k-1)J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)J(K,I+1)=r; %DF/U(i,j,k+1)J(K,I+N)=r; %DF/U(i,j+1,k)J(K,I+N*N)=r; %DF/U(i+1,j,k)endif i==1 && j==1 && k~=1 && k==NJ(K,I-1)=q+r; %DF/U(i,j,k-1)J(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)J(K,I+N)=r; %DF/U(i,j+1,k)J(K,I+N*N)=r; %DF/U(i+1,j,k)end%i==1 && j~=1 && k==1 && j~=Nif i==1 && j~=1 && k==1 && j~=NJ(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)J(K,I+1)=r; %DF/U(i,j,k+1)J(K,I+N)=r; %DF/U(i,j+1,k)J(K,I-N)=q+r; %DF/U(i,j-1,k)J(K,I+N*N)=r; %DF/U(i+1,j,k)endif i==1 && j~=1 && k==1 && j==NJ(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)J(K,I+1)=r; %DF/U(i,j,k+1)J(K,I-N)=q+r; %DF/U(i,j-1,k)J(K,I+N*N)=r; %DF/U(i+1,j,k)end%i~=1 && j==1 && k==if i~=1 && j==1 && k==1 && i~=NJ(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)J(K,I+1)=r; %DF/U(i,j,k+1)J(K,I+N)=r; %DF/U(i,j+1,k)J(K,I+N*N)=r; %DF/U(i+1,j,k)J(K,I-N*N)=q+r; %DF/U(i-1,j,k)endif i~=1 && j==1 && k==1 && i==NJ(K,I)=-3*q-(6*r+1); %DF/U(i,j,k)J(K,I+1)=r; %DF/U(i,j,k+1)J(K,I+N)=r; %DF/U(i,j+1,k)J(K,I-N*N)=q+r; %DF/U(i-1,j,k)endend% Solving the system of equations J*dU=-F, we can use LU decomposition to save% computational resources in robust systemdU=(J(-F));EPS=max(abs(dU(:,1)));a=0;for i=1:N %xfor j=1:N %yfor k=1:N %za=a+1;if i~=1 && j~=1 && k~=1 && i~=N && j~=N && k~=NU(i,j,k,n+1)=U(i,j,k,n+1)+dU(a);endendendenditeration=iteration+1;endend %% Plots and subrotines for 3D plottingm = createMesh3D(N-2,N-2,N-2,Lx,Ly,Lz); % create the meshNumerical=CellVariable(m, U(:,:,:,M));Analytical=CellVariable(m, UA(:,:,:,M));Error=CellVariable(m,abs(U(:,:,:,M)-UA(:,:,:,M)));figurevisualizeCells3D(Numerical);title('Numerical solution at t=0.25')figurevisualizeCells3D(Analytical);title('Analytical solution at t=0.25')figurevisualizeCells3D(Error);title('Absolute Error at t=0.25')function visualizeCells3D(phi)% Modified from 2012-2016 Ali Akbar Eftekhariphi.value = phi.value(2:end-1,2:end-1,2:end-1);[X,Y,Z]=meshgrid(phi.domain.cellcenters.y, phi.domain.cellcenters.x, ...phi.domain.cellcenters.z);phi.value(1)=phi.value(1)+eps; % to avoid an strange error for assigning color limitsSx = [phi.domain.cellcenters.x(1) phi.domain.cellcenters.x(end)];Sy = [phi.domain.cellcenters.y(1) phi.domain.cellcenters.y(end)];Sz = [phi.domain.cellcenters.z(1) phi.domain.cellcenters.z(end)];slice(X,Y,Z, phi.value, Sy, Sx, Sz);xlabel('[y vlaues]'); % this is correct [matrix not rotated]ylabel('[x vlaues]'); % this is correct [matrix not rotated]zlabel('[z vlaues]');axis equal tightcolorbarend function MS = createMesh3D(varargin)% Modified from Ali Akbar EftekhariNx=varargin{1};Ny=varargin{2};Nz=varargin{3};Width=varargin{4};Height=varargin{5};Depth=varargin{6};% cell size is dxdx =0.125;dy = 0.125;dz = 0.125;G=reshape(1:(Nx+2)*(Ny+2)*(Nz+2), Nx+2, Ny+2, Nz+2);CellSize.x= dx*ones(Nx+2,1);CellSize.y= dy*ones(Ny+2,1);CellSize.z= dz*ones(Nz+2,1);CellLocation.x= [1:Nx]'*dx;CellLocation.y= [1:Ny]'*dy;CellLocation.z= [1:Nz]'*dz;FaceLocation.x= [0:Nx]'*dx;FaceLocation.y= [0:Ny]'*dy;FaceLocation.z= [0:Nz]'*dz;c=G([1,end], [1,end], [1, end]);e1=G([1, end], [1, end], 2:Nz+1);e2=G([1, end], 2:Ny+1, [1, end]);e3=G(2:Nx+1, [1, end], [1, end]);MS=createMesh3D(3,[Nx,Ny, Nz],CellSize,CellLocation,FaceLocation, c(:),[e1(:); e2(:); e3(:)]);end % Used for 3D plottingclassdef HW7% From 2012-2016 Ali Akbar Eftekharipropertiesdimensiondimscellsizecellcentersfacecenterscornersedgesendmethodsfunction meshVar = HW7(dimension, dims, cellsize, ...cellcenters, facecenters, corners, edges)if nargin>0meshVar.dimension = dimension;meshVar.dims = dims;meshVar.cellsize = cellsize;meshVar.cellcenters = cellcenters;meshVar.facecenters = facecenters;meshVar.corners= corners;meshVar.edges= edges;endendendend Computer ScienceEngineering & TechnologyJava ProgrammingMATH 216Share Question
