How can we make the following spectral relaxation code to run?
Show older comments
The dimensionless sisko models to be solved by the spectral relaxation method are:
function [D, x] = cheb(N)
if N == 0
D = 0;
x = 1;
return
end
x = cos(pi * (0:N) / N)';
c = [2; ones(N - 1, 1); 2] .* (-1).^(0:N)';
X = repmat(x, 1, N + 1);
dX = X - X';
D = (c * (1 ./ c)')./(dX + (eye(N + 1)));
D = D - diag(sum(D, 2));
end
clear; %clc;
N = 120;
[D,x] = cheb(N);
% Parameter initialization
% A = 1; %Material parameter
M = 1; %Magnetic field
n = 1; % Power law index
m = 1; %Velocity power
gamma = 1; %Chemical reaction parameter
Nt = 0.1; %Thermophoresis number
Nb = 0.1; %Brownian motion
GrT = 1; %Thermal Grashof number
GrC = 2; %Concentration Grashof number
Rd = 1; %Radiation parameter
Sc = 1; %Schmidt number
Pr = 3; %Prandtl number
Ec = 1; %Eckert number
%Sk=linspace(0,0.5,6); % Domain for the Skin friction coeeficient,Nusselt number and sherwood number
% Constant calculation
Costant = (m*(2*n - 1) + 1) / (n + 1);
iterations = 30;
tol = 1e-8;
tic
p = 1; %counter
Vec =[ 1 2 3 4];
for mm = 1: length(Vec)
A= Vec(mm);
Lmin = 1; MaxL = 10;
L = Lmin;
while L < MaxL
scale = 2/L;
D1 = scale*D; D2 = (scale^2*D^2) ; D3 = (scale^3*D^3);
y = L*(x+1)/2; I = eye(N+1);
gr = exp(-y);hr = exp(-y); thetar = exp(-y); phir = exp(-y);
fr = 1- exp(-y);
it = 1;
thetanorm = 1; phinorm = 1;hnorm = 1; gnorm = 1;
NORMS = [gnorm hnorm phinorm thetanorm];
while max(NORMS) > tol && it < iterations
%=================================================================================
%Momentum equation
A1 = (A+n*(-D1*diag(gr)).^(n-1))*D2+Costant*fr*D1-M.*I; %Left hand side A1*g_{r+1}=B1
A1(N+1,:) = 0; A1(N+1,N+1) = 1;
A1(1,:) = 0; A1(1,1) = 1;
B1 =m*gr^2 -GrT*thetar - GrC*phir; %Right hand side
B1(1) = 0; B1(N+1) = 1;
grt = A1\B1;
gnorm = norm(grt-gr,inf);
gr = grt;
gr1 = D1*gr;
%================================================================================
%Energy equation
A2 =(1/Pr)*(1+(4/3)*Rd)*D2+Costant*fr*D1+Nb*(D1*phir)*D1;
A2(N+1,:) = 0; A2(N+1,N+1) = 1;
A2(1,:) = 0; A2(1,1) = 1;
B1 = -M*Ec*(gr)^2-Nt*(D1*thetar)^2-Ec*A*(D1*gr)^(2)-Ec*(-D1*gr)^(n+1);
B1(1) = 0; B1(N+1) = 1;
thetart = A2\B1;
hnorm = norm(thetart-thetar,inf);
thetar = thetart;
thetar1 = D1*thetar;
%=======================================================================================
%Concentration equation
A3 = D2+Sc*Costant*fr*D-gamma*I ;
A3(N+1,:) = 0; A3(N+1,N+1) = 1;
A3(1,:) = 0; A3(1,1) = 1;
B3 = -(Nt/Nb)*(D2*thetar);
B3(1) = 0;
B3(N+1) = 1;
phirt = A3\B3;
thetanorm = norm(phirt-phir,inf);
phir = phirt;
phir1 = D1*phir;
%--------------------------------------------------------------------------
% A4 = D2 + Le*diag(fr)*D1;
% A4(N+1,:) = 0; A4(N+1,N+1) = 1;
% A4(1,:) = 0; A4(1,1) = 1;
% B4 = -(Nt/Nb)*(D2*thetar);
% B4(1) = 0;
% B4(N+1) = 1;
% phirt = inv(A4)*B4;
% phinorm = norm(phirt-phir,inf);
% phir = phirt;
% phir1 = D1*phir;
%--------------------------------------------------------------------------
fprintf('%10.0f\t %10.10f\t %10.10f\t %10.10f\n', it,-gr1(N+1),-thetar1(N+1),-phir1(N+1));
it = it + 1;
NORMS = [gnorm hnorm phinorm thetanorm];
Edecoupled(it-1) = max(NORMS);
end
thetarr(:,p+1) = thetar1(N+1); %f'(eta) after convergence
phirr(:,p+1) = phir1(N+1);
grr(:,p+1) = gr1(N+1);
dtheta = abs(thetarr(:,p+1) - thetarr(:,p));
dphi = abs(phirr(:,p+1)-phirr(:,p));
dg = abs(grr(:,p+1)-grr(:,p));
dL = [dtheta dphi dg];
if all(dL) < tol
% if dtheta <tol
break
end
L = L + 1;
p = p + 1; %Counter
end
switch mm
case 1
figure(1)
plot(1:length(Edecoupled),log(Edecoupled),'ko-','markersize',5,'MarkerFaceColor','k')
clear Edecoupled;
hold on
figure(2)
plot(y,gr,'k','linewidth',1.5);
hold on
figure(3)
plot(y,thetar,'k','linewidth',1.5);
hold on
figure(4)
plot(y,phir,'k','linewidth',1.5);
hold on
case 2
figure(1)
plot(1:length(Edecoupled),log(Edecoupled),'k*-','markersize',5,'MarkerFaceColor','k')
clear Edecoupled;
figure(2)
plot(y,gr,'r:','linewidth',1.5);
hold on
figure(3)
plot(y,thetar,'r:','linewidth',1.5);
hold on
figure(4)
plot(y,phir,'r:','linewidth',1.5);
hold on
case 3
figure(1)
plot(1:length(Edecoupled),log(Edecoupled),'k^-','markersize',5,'MarkerFaceColor','k')
clear Edecoupled;
figure(2)
plot(y,gr,'g--','linewidth',1.5);
hold on
figure(3)
plot(y,thetar,'g--','linewidth',1.5);
hold on
figure(4)
plot(y,phir,'g--','linewidth',1.5);
hold on
case 4
figure(1)
plot(1:length(Edecoupled),log(Edecoupled),'kd-','markersize',5,'MarkerFaceColor','k')
clear Edecoupled;
figure(2)
plot(y,gr,'b--','linewidth',1.5);
hold on
figure(3)
plot(y,thetar,'b--','linewidth',1.5);
hold on
figure(4)
plot(y,phir,'b--','linewidth',1.5);
hold on
case 5
figure(1)
plot(1:length(Edecoupled),log(Edecoupled),'kp-','markersize',5,'MarkerFaceColor','k')
clear Edecoupled;
figure(2)
plot(y,gr,'b---','linewidth',1.5);
hold on
figure(3)
plot(y,thetar,'b---','linewidth',1.5);
hold on
figure(4)
plot(y,phir,'b---','linewidth',1.5);
draw
hold on
% case 5
% figure(1)
%plot(1:length(Edecoupled),log(Edecoupled),'kp-','markersize',5,'MarkerFaceColor','k')
% clear Edecoupled;
end
figure(1)
xlabel('Iterations','FontSize',12)
ylabel('$\ln(E_d)$','Interpreter','Latex','FontSize',14)
legend([repmat('M=',length(MM),1) num2str(MM')])
figure(2)
legend([repmat('M= ',length(MM),1) num2str(MM')])
xlabel('$\eta$','Interpreter','Latex','FontSize',14)
ylabel('$f''(\eta)$','Interpreter','Latex','Rotation',90,'FontSize',14)
axis([0 10 0 1])
figure(3)
legend([repmat('M= ',length(MM),1) num2str(MM')])
xlabel('$\eta$','Interpreter','Latex','FontSize',14)
ylabel('$\theta(\eta)$','Interpreter','Latex','Rotation',90,'FontSize',14)
axis([0 10 0 1])
figure(4)
legend([repmat('M= ',length(MM),1) num2str(MM')])
xlabel('$\eta$','Interpreter','Latex','FontSize',14)
ylabel('$\phi(\eta)$','Interpreter','Latex','Rotation',90,'FontSize',14)
axis([0 10 0 1])
%fprintf('%2.0f\t %2.0f\t %2.0f\t %2.0f\t %2.0f\t %2.0f\t %2.0f\t %10.8f\t %10.8f\n',G_r,K_r,it,L,N,pr1(N+1),thetar1(N+1));
end
toc
Accepted Answer
Walter Roberson
on 22 Apr 2026
The actual error message is
Error using *
Incorrect dimensions for matrix multiplication. Check that the number of columns in the first matrix matches the number of rows in the second matrix. To operate on each element
of the matrix individually, use TIMES (.*) for elementwise multiplication.
Error in
sisko (line 57)
A1 = (A+n*(-D1*diag(gr)).^(n-1))*D2+Costant*fr*D1-M.*I; %Left hand side A1*g_{r+1}=B1
(A+n*(-D1*diag(gr)).^(n-1))*D2 is 121 x 121
Costant is 1 x 1
fr is 121 x 1
D1 is 121 x 121
M.*I is 121 x 121
So you have the problem that you are using the * operation between a 121 x 1 and a 121 x 121.
Now,
fr = 1- exp(-y);
and
y = L*(x+1)/2;
where L is scalar and x is 121 x 1, so y is 121 x 1.
In turn
[D,x] = cheb(N);
Inside cheb()
x = cos(pi * (0:N) / N)';
which clearly forms a row vector and then transposes it into a column vector. So x is going to be a column vector, not an array, so it makes sense that by the time fr is formed that fr is a column vector, and thus that it is an error to use the * operation between 121 x 1 column vector fr and 121 x 121 array D1.
You need to decide if you want is fr'*D1 -- which would make it a 1 x 121 that is * by 121 x 121, giving a 1 x 121 result. It is permitted to + between the 121 x 121 array before that on the A1 line, but it is not necessarily what you want.
You need to decide if what you want is fr.*D1 -- which would make it 121 x 1 .* 121 x 121 which would produce a 121 x 121 result.
But possibly you are expecting fr to be 121 x 121 at that point. If so then you are not constructing it properly.
3 Comments
But is there any way we can imrpove the code:
clear; %clc;
N = 120;
[D,x] = cheb(N);
% Parameter initialization
% A = 1; %Material parameter
M = 1; %Magnetic field
n = 1; % Power law index
m = 1; %Velocity power
gamma = 1; %Chemical reaction parameter
Nt = 0.1; %Thermophoresis number
Nb = 0.1; %Brownian motion
GrT = 1; %Thermal Grashof number
GrC = 2; %Concentration Grashof number
Rd = 1; %Radiation parameter
Sc = 1; %Schmidt number
Pr = 3; %Prandtl number
Ec = 1; %Eckert number
%Sk=linspace(0,0.5,6); % Domain for the Skin friction coeeficient,Nusselt number and sherwood number
% Constant calculation
Costant = (m*(2*n - 1) + 1) / (n + 1);
iterations = 30;
tol = 1e-8;
tic
p = 1; %counter
Vec =[ 1 2 3 4];
for mm = 1: length(Vec)
A= Vec(mm);
Lmin = 1; MaxL = 10;
L = Lmin;
while L < MaxL
scale = 2/L;
D1 = scale*D; D2 = (scale^2*D^2) ; D3 = (scale^3*D^3);
y = L*(x+1)/2; I = eye(N+1);
gr = exp(-y);hr = exp(-y); thetar = exp(-y); phir = exp(-y);
fr = 1- exp(-y);
it = 1;
thetanorm = 1; phinorm = 1;hnorm = 1; gnorm = 1;
NORMS = [gnorm hnorm phinorm thetanorm];
while max(NORMS) > tol && it < iterations
%=================================================================================
%Momentum equation
A1 = (A+n*(-D1*diag(gr)).^(n-1))*D2+Costant*fr.*D1-M.*I; %Left hand side A1*g_{r+1}=B1
A1(N+1,:) = 0; A1(N+1,N+1) = 1;
A1(1,:) = 0; A1(1,1) = 1;
B1 =m*gr.^2 -GrT*thetar - GrC*phir; %Right hand side
B1(1) = 0; B1(N+1) = 1;
grt = A1\B1;
gnorm = norm(grt-gr,inf);
gr = grt;
gr1 = D1*gr;
%================================================================================
%Energy equation
A2 =(1/Pr)*(1+(4/3)*Rd)*D2+Costant*fr.*D1+Nb*(D1*phir).*D1;
A2(N+1,:) = 0; A2(N+1,N+1) = 1;
A2(1,:) = 0; A2(1,1) = 1;
B2 = -M*Ec*(gr).^2-Nt*(D1*thetar).^2-Ec*A*(D1*gr).^(2)-Ec*(-D1*gr).^(n+1);
B2(1) = 0; B2(N+1) = 1;
thetart = A2\B2;
hnorm = norm(thetart-thetar,inf);
thetar = thetart;
thetar1 = D1*thetar;
%=======================================================================================
%Concentration equation
A3 = D2+Sc*Costant*fr.*D-gamma*I ;
A3(N+1,:) = 0; A3(N+1,N+1) = 1;
A3(1,:) = 0; A3(1,1) = 1;
B3 = -(Nt/Nb)*(D2*thetar);
B3(1) = 0;
B3(N+1) = 1;
phirt = A3\B3;
thetanorm = norm(phirt-phir,inf);
phir = phirt;
phir1 = D1*phir;
%--------------------------------------------------------------------------
% A4 = D2 + Le*diag(fr)*D1;
% A4(N+1,:) = 0; A4(N+1,N+1) = 1;
% A4(1,:) = 0; A4(1,1) = 1;
% B4 = -(Nt/Nb)*(D2*thetar);
% B4(1) = 0;
% B4(N+1) = 1;
% phirt = inv(A4)*B4;
% phinorm = norm(phirt-phir,inf);
% phir = phirt;
% phir1 = D1*phir;
%--------------------------------------------------------------------------
fprintf('%10.0f\t %10.10f\t %10.10f\t %10.10f\n', it,-gr1(N+1),-thetar1(N+1),-phir1(N+1));
it = it + 1;
NORMS = [gnorm hnorm phinorm thetanorm];
Edecoupled(it-1) = max(NORMS);
end
thetarr(:,p+1) = thetar1(N+1); %f'(eta) after convergence
phirr(:,p+1) = phir1(N+1);
grr(:,p+1) = gr1(N+1);
dtheta = abs(thetarr(:,p+1) - thetarr(:,p));
dphi = abs(phirr(:,p+1)-phirr(:,p));
dg = abs(grr(:,p+1)-grr(:,p));
dL = [dtheta dphi dg];
if all(dL) < tol
% if dtheta <tol
break
end
L = L + 1;
p = p + 1; %Counter
end
switch mm
case 1
figure(1)
plot(1:length(Edecoupled),log(Edecoupled),'ko-','markersize',5,'MarkerFaceColor','k')
clear Edecoupled;
hold on
figure(2)
plot(y,gr,'k','linewidth',1.5);
hold on
figure(3)
plot(y,thetar,'k','linewidth',1.5);
hold on
figure(4)
plot(y,phir,'k','linewidth',1.5);
hold on
case 2
figure(1)
plot(1:length(Edecoupled),log(Edecoupled),'k*-','markersize',5,'MarkerFaceColor','k')
clear Edecoupled;
figure(2)
plot(y,gr,'r:','linewidth',1.5);
hold on
figure(3)
plot(y,thetar,'r:','linewidth',1.5);
hold on
figure(4)
plot(y,phir,'r:','linewidth',1.5);
hold on
case 3
figure(1)
plot(1:length(Edecoupled),log(Edecoupled),'k^-','markersize',5,'MarkerFaceColor','k')
clear Edecoupled;
figure(2)
plot(y,gr,'g--','linewidth',1.5);
hold on
figure(3)
plot(y,thetar,'g--','linewidth',1.5);
hold on
figure(4)
plot(y,phir,'g--','linewidth',1.5);
hold on
case 4
figure(1)
plot(1:length(Edecoupled),log(Edecoupled),'kd-','markersize',5,'MarkerFaceColor','k')
clear Edecoupled;
figure(2)
plot(y,gr,'b--','linewidth',1.5);
hold on
figure(3)
plot(y,thetar,'b--','linewidth',1.5);
hold on
figure(4)
plot(y,phir,'b--','linewidth',1.5);
hold on
case 5
figure(1)
plot(1:length(Edecoupled),log(Edecoupled),'kp-','markersize',5,'MarkerFaceColor','k')
clear Edecoupled;
figure(2)
plot(y,gr,'b---','linewidth',1.5);
hold on
figure(3)
plot(y,thetar,'b---','linewidth',1.5);
hold on
figure(4)
plot(y,phir,'b---','linewidth',1.5);
draw
hold on
% case 5
% figure(1)
%plot(1:length(Edecoupled),log(Edecoupled),'kp-','markersize',5,'MarkerFaceColor','k')
% clear Edecoupled;
end
figure(1)
xlabel('Iterations','FontSize',12)
ylabel('$\ln(E_d)$','Interpreter','Latex','FontSize',14)
legend([repmat('A=',length(Vec),1) num2str(Vec')])
figure(2)
legend([repmat('A= ',length(Vec),1) num2str(Vec')])
xlabel('$\eta$','Interpreter','Latex','FontSize',14)
ylabel('$f''(\eta)$','Interpreter','Latex','Rotation',90,'FontSize',14)
axis([0 10 0 1])
figure(3)
legend([repmat('A= ',length(Vec),1) num2str(Vec')])
xlabel('$\eta$','Interpreter','Latex','FontSize',14)
ylabel('$\theta(\eta)$','Interpreter','Latex','Rotation',90,'FontSize',14)
axis([0 10 0 1])
figure(4)
legend([repmat('A= ',length(Vec),1) num2str(Vec')])
xlabel('$\eta$','Interpreter','Latex','FontSize',14)
ylabel('$\phi(\eta)$','Interpreter','Latex','Rotation',90,'FontSize',14)
axis([0 10 0 1])
%fprintf('%2.0f\t %2.0f\t %2.0f\t %2.0f\t %2.0f\t %2.0f\t %2.0f\t %10.8f\t %10.8f\n',G_r,K_r,it,L,N,pr1(N+1),thetar1(N+1));
end
toc
function [D, x] = cheb(N)
if N == 0
D = 0;
x = 1;
return
end
x = cos(pi * (0:N) / N)';
c = [2; ones(N - 1, 1); 2] .* (-1).^(0:N)';
X = repmat(x, 1, N + 1);
dX = X - X';
D = (c * (1 ./ c)')./(dX + (eye(N + 1)));
D = D - diag(sum(D, 2));
end
so that it may plot and converge?
If we are confirming the graphs they must look like the ones on the bvp4c:
sisko9()
function sisko9
% Parameter initialization
% A = 1; %Material parameter
M = 1; %Magnetic field
n = 1; % Power law index
m = 1; %Velocity power
gamma = 1; %Chemical reaction parameter
Nt = 0.1; %Thermophoresis number
Nb = 0.1; %Brownian motion
GrT = 1; %Thermal Grashof number
GrC = 2; %Concentration Grashof number
Rd = 1; %Radiation parameter
Sc = 1; %Schmidt number
Pr = 3; %Prandtl number
Ec = 1; %Eckert number
%Sk=linspace(0,0.5,6); % Domain for the Skin friction coeeficient,Nusselt number and sherwood number
% Constant calculation
Costant = (m*(2*n - 1) + 1) / (n + 1);
% System of ODE
function dydx = odefun(~, y)
%f = y(1), f' = y(2), f'' = y(3), theta = y(4), theta' = y(5),
%phi=y(6), phi'=y(7)
%= y(6), phi' = y(7), phi''=y(8)
dydx = zeros(7, 1);
dydx(1) = y(2);
dydx(2) = y(3);
dydx(3) = (-Costant*y(1)*y(3) + m*y(2)^2 - GrT*y(4) - GrC*y(6) + M*y(2))./(A(i) + n*(-y(3))^(n-1));
dydx(4) = y(5);
dydx(5) = (-Costant*y(1)*y(5) - Nb*y(5)*y(7) - Nt*y(5)^2 - Ec.*A(i).*y(3).^2 - Ec*(-y(3))^(n+1) - M*Ec*y(2)^2)./((1+4/3*Rd)/Pr);
dydx(6) = y(7);
dydx(7) = -Sc*Costant*y(1)*y(7) - (Nt/Nb)*dydx(5) + gamma*y(6);
end
% Nested function: Boundary conditions
function res = bcfun(ya, yb)
res = zeros(7, 1);
res(1) = ya(2)-1; % y(2) = 1 at eta = 0
res(2) = ya(1); % y(1) = 0 at eta = 0
res(3) = ya(4)-1; % y(4) = 1 at eta = 0
res(4) = Nb*ya(7) + Nt*ya(5); % Boundary condition at eta = 0
res(5) = yb(2); % y(2) = 0 at eta = 10
res(6) = yb(4); % y(4) = 0 at eta = 10
res(7) = yb(6); % y(6) = 0 at eta = 10
end
% Classical variables
%
%==============================================================================
solinit=bvpinit(linspace(0,10,10),[0 0 0 0 0 0 0]);
A=[1 2 3 4 ];
C={'r','g','k','b'}; % Original models colors for parameter variation
%Cl={'r','g','b','k'}; % Skin friction, Nusselt number and Sherwood number colors
for i=1:4
% options=bvpset('RelTol',1e-10,'AbsTol',1e-10);
sol=bvp4c(@odefun,@bcfun,solinit);
eta=linspace(0,10,400);
y=deval(sol,eta);
Cf_x=A(i)*y(3,1)-(-y(3,1)); % Skin friction coefficient
% Nu_x=-(1+(4/3)*Rd)*y(5,1); %Nusselt number
% Sh_x=-y(7,1); %Sherwood number
%=====================================================================================
figure(1);
plot(eta,y(1,:),C{i},Linestyle='-',LineWidth=2);
xlabel("\eta");
ylabel("f(\eta)");
% legend('A=1','A=2','A=3','A=4')
hold on
%=====================================================================================
figure(2);
plot(eta,y(2,:),C{i},LineStyle="-",LineWidth=2);
xlabel("\eta");
ylabel("f'(\eta)");
%legend('A=1','A=2','A=3','A=4');
hold on
%=====================================================================================
figure(3);
plot(eta,y(3,:),C{i},LineStyle="-",LineWidth=2);
xlabel("\eta");
ylabel("f''(\eta)");
%legend('A=1','A=2','A=3','A=4');
hold on
%====================================================================================
figure(4);
plot(eta,y(4,:),C{i},LineStyle="-",LineWidth=2);
xlabel("\eta");
ylabel("\theta(\eta)");
%legend('A=1','A=2','A=3','A=4');
hold on
%====================================================================================
figure(5);
plot(eta,y(5,:),C{i},LineStyle="-",LineWidth=2);
xlabel("\eta");
ylabel("\theta'(\eta)");
%legend('A=1','A=2','A=3','A=4');
hold on
%=====================================================================
figure(6);
plot(eta,y(6,:),C{i},LineStyle="-",LineWidth=2);
xlabel("\eta");
ylabel("\phi(\eta)");
%legend('A=1','A=2','A=3','A=4');
hold on
%=========================================================================
figure(7);
plot(eta,y(7,:),C{i},LineStyle="-",LineWidth=2);
xlabel("\eta");
ylabel("\phi''(\eta)");
%legend('A=1','A=2','A=3','A=4');
hold on
figure(8);
plot(A(i),Cf_x,strcat(C{i},'*'));
xlabel("A");
ylabel('$(1/2)Re_{b}^{-1/(n+1)}Cf_{x}$','Interpreter','latex')
hold on
% %============================================================
% figure(9);
% plot(eta,Nu_x,C{i},LineStyle='-',LineWidth=2);
% xlabel("\eta");
% ylabel('$Re_{b}^{-1/(n+1)}Nu_{x}$');
%
% hold on
% figure(10);
% plot(eta,Sh_x,C{i},LineStyle="-",LineWidth=2);
% xlabel("\eta");
% ylabel('$Re_{b}^{-1/(n+1)}Sh_{x}$');
%hold on
end
figure(1);
legend([repmat('A= ',length(A),1) num2str(A')]);
figure(2);
legend([repmat('A= ',length(A),1) num2str(A')]);
figure(3);
legend([repmat('A= ',length(A),1) num2str(A')]);
figure(4);
legend([repmat('A= ',length(A),1) num2str(A')]);
figure(5);
legend([repmat('A= ',length(A),1) num2str(A')]);
figure(6);
legend([repmat('A= ',length(A),1) num2str(A')]);
figure(7);
legend([repmat('A= ',length(A),1) num2str(A')]);
figure(8);
legend([repmat('A= ',length(A),1) num2str(A')]);
% legend('A=1','A=2','A=3','A=4');
% figure(9);
% legend('A=1','A=2','A=3','A=4');
% axis([0 1 0 1]);
% figure(10);
% legend('A=1','A=2','A=3','A=4');
% for i=
% %fprintf('%.10f\n',Cf_x);
% end
%end
end
More Answers (2)
I modified A1, B1, A2, B1 and A3 to get compatible dimensions for matrix multiplications without knowing whether these modifications are correct from the mathematical point of view.
The code now works technically, but doesn't converge.
Since you don't explain the mathematical method you try to implement (and I've not yet heard of it), I cannot help further.
function [D, x] = cheb(N)
if N == 0
D = 0;
x = 1;
return
end
x = cos(pi * (0:N) / N)';
c = [2; ones(N - 1, 1); 2] .* (-1).^(0:N)';
X = repmat(x, 1, N + 1);
dX = X - X';
D = (c * (1 ./ c)')./(dX + (eye(N + 1)));
D = D - diag(sum(D, 2));
end
clear; %clc;
N = 120;
[D,x] = cheb(N);
% Parameter initialization
% A = 1; %Material parameter
M = 1; %Magnetic field
n = 1; % Power law index
m = 1; %Velocity power
gamma = 1; %Chemical reaction parameter
Nt = 0.1; %Thermophoresis number
Nb = 0.1; %Brownian motion
GrT = 1; %Thermal Grashof number
GrC = 2; %Concentration Grashof number
Rd = 1; %Radiation parameter
Sc = 1; %Schmidt number
Pr = 3; %Prandtl number
Ec = 1; %Eckert number
%Sk=linspace(0,0.5,6); % Domain for the Skin friction coeeficient,Nusselt number and sherwood number
% Constant calculation
Costant = (m*(2*n - 1) + 1) / (n + 1);
iterations = 30;
tol = 1e-8;
tic
p = 1; %counter
Vec =[ 1 2 3 4];
for mm = 1: length(Vec)
A= Vec(mm);
Lmin = 1; MaxL = 10;
L = Lmin;
while L < MaxL
scale = 2/L;
D1 = scale*D; D2 = (scale^2*D^2) ; D3 = (scale^3*D^3);
y = L*(x+1)/2; I = eye(N+1);
gr = exp(-y);hr = exp(-y); thetar = exp(-y); phir = exp(-y);
fr = 1- exp(-y);
it = 1;
thetanorm = 1; phinorm = 1;hnorm = 1; gnorm = 1;
NORMS = [gnorm hnorm phinorm thetanorm];
while max(NORMS) > tol && it < iterations
%=================================================================================
%Momentum equation
% A1 = (A+n*(-D1*diag(gr)).^(n-1))*D2+Costant*fr*D1-M.*I; %Left hand side A1*g_{r+1}=B1
A1 = (A+n*(-D1*diag(gr)).^(n-1))*D2+Costant*D1*fr-M.*I; %Left hand side A1*g_{r+1}=B1
A1(N+1,:) = 0; A1(N+1,N+1) = 1;
A1(1,:) = 0; A1(1,1) = 1;
% B1 =m*gr^2 -GrT*thetar - GrC*phir; %Right hand side
B1 =m*gr.^2 -GrT*thetar - GrC*phir; %Right hand side
B1(1) = 0; B1(N+1) = 1;
grt = A1\B1;
gnorm = norm(grt-gr,inf);
gr = grt;
gr1 = D1*gr;
%================================================================================
%Energy equation
% A2 =(1/Pr)*(1+(4/3)*Rd)*D2+Costant*fr*D1+Nb*(D1*phir)*D1;
A2 =(1/Pr)*(1+(4/3)*Rd)*D2+Costant*D1*fr+Nb*D1*(D1*phir);
A2(N+1,:) = 0; A2(N+1,N+1) = 1;
A2(1,:) = 0; A2(1,1) = 1;
% B1 = -M*Ec*(gr)^2-Nt*(D1*thetar)^2-Ec*A*(D1*gr)^(2)-Ec*(-D1*gr)^(n+1);
B1 = -M*Ec*(gr).^2-Nt*(D1*thetar).^2-Ec*A*(D1*gr).^(2)-Ec*(-D1*gr).^(n+1);
B1(1) = 0; B1(N+1) = 1;
thetart = A2\B1;
hnorm = norm(thetart-thetar,inf);
thetar = thetart;
thetar1 = D1*thetar;
%=======================================================================================
%Concentration equation
% A3 = D2+Sc*Costant*fr*D-gamma*I ;
A3 = D2+Sc*Costant*D*fr-gamma*I;
A3(N+1,:) = 0; A3(N+1,N+1) = 1;
A3(1,:) = 0; A3(1,1) = 1;
B3 = -(Nt/Nb)*(D2*thetar);
B3(1) = 0;
B3(N+1) = 1;
phirt = A3\B3;
thetanorm = norm(phirt-phir,inf);
phir = phirt;
phir1 = D1*phir;
%--------------------------------------------------------------------------
% A4 = D2 + Le*diag(fr)*D1;
% A4(N+1,:) = 0; A4(N+1,N+1) = 1;
% A4(1,:) = 0; A4(1,1) = 1;
% B4 = -(Nt/Nb)*(D2*thetar);
% B4(1) = 0;
% B4(N+1) = 1;
% phirt = inv(A4)*B4;
% phinorm = norm(phirt-phir,inf);
% phir = phirt;
% phir1 = D1*phir;
%--------------------------------------------------------------------------
fprintf('%10.0f\t %10.10f\t %10.10f\t %10.10f\n', it,-gr1(N+1),-thetar1(N+1),-phir1(N+1));
it = it + 1;
NORMS = [gnorm hnorm phinorm thetanorm];
Edecoupled(it-1) = max(NORMS);
end
thetarr(:,p+1) = thetar1(N+1); %f'(eta) after convergence
phirr(:,p+1) = phir1(N+1);
grr(:,p+1) = gr1(N+1);
dtheta = abs(thetarr(:,p+1) - thetarr(:,p));
dphi = abs(phirr(:,p+1)-phirr(:,p));
dg = abs(grr(:,p+1)-grr(:,p));
dL = [dtheta dphi dg];
if all(dL) < tol
% if dtheta <tol
break
end
L = L + 1;
p = p + 1; %Counter
end
switch mm
case 1
figure(1)
plot(1:length(Edecoupled),log(Edecoupled),'ko-','markersize',5,'MarkerFaceColor','k')
clear Edecoupled;
hold on
figure(2)
plot(y,gr,'k','linewidth',1.5);
hold on
figure(3)
plot(y,thetar,'k','linewidth',1.5);
hold on
figure(4)
plot(y,phir,'k','linewidth',1.5);
hold on
case 2
figure(1)
plot(1:length(Edecoupled),log(Edecoupled),'k*-','markersize',5,'MarkerFaceColor','k')
clear Edecoupled;
figure(2)
plot(y,gr,'r:','linewidth',1.5);
hold on
figure(3)
plot(y,thetar,'r:','linewidth',1.5);
hold on
figure(4)
plot(y,phir,'r:','linewidth',1.5);
hold on
case 3
figure(1)
plot(1:length(Edecoupled),log(Edecoupled),'k^-','markersize',5,'MarkerFaceColor','k')
clear Edecoupled;
figure(2)
plot(y,gr,'g--','linewidth',1.5);
hold on
figure(3)
plot(y,thetar,'g--','linewidth',1.5);
hold on
figure(4)
plot(y,phir,'g--','linewidth',1.5);
hold on
case 4
figure(1)
plot(1:length(Edecoupled),log(Edecoupled),'kd-','markersize',5,'MarkerFaceColor','k')
clear Edecoupled;
figure(2)
plot(y,gr,'b--','linewidth',1.5);
hold on
figure(3)
plot(y,thetar,'b--','linewidth',1.5);
hold on
figure(4)
plot(y,phir,'b--','linewidth',1.5);
hold on
case 5
figure(1)
plot(1:length(Edecoupled),log(Edecoupled),'kp-','markersize',5,'MarkerFaceColor','k')
clear Edecoupled;
figure(2)
plot(y,gr,'b---','linewidth',1.5);
hold on
figure(3)
plot(y,thetar,'b---','linewidth',1.5);
hold on
figure(4)
plot(y,phir,'b---','linewidth',1.5);
draw
hold on
% case 5
% figure(1)
%plot(1:length(Edecoupled),log(Edecoupled),'kp-','markersize',5,'MarkerFaceColor','k')
% clear Edecoupled;
end
figure(1)
xlabel('Iterations','FontSize',12)
ylabel('$\ln(E_d)$','Interpreter','Latex','FontSize',14)
legend([repmat('M=',length(MM),1) num2str(MM')])
figure(2)
legend([repmat('M= ',length(MM),1) num2str(MM')])
xlabel('$\eta$','Interpreter','Latex','FontSize',14)
ylabel('$f''(\eta)$','Interpreter','Latex','Rotation',90,'FontSize',14)
axis([0 10 0 1])
figure(3)
legend([repmat('M= ',length(MM),1) num2str(MM')])
xlabel('$\eta$','Interpreter','Latex','FontSize',14)
ylabel('$\theta(\eta)$','Interpreter','Latex','Rotation',90,'FontSize',14)
axis([0 10 0 1])
figure(4)
legend([repmat('M= ',length(MM),1) num2str(MM')])
xlabel('$\eta$','Interpreter','Latex','FontSize',14)
ylabel('$\phi(\eta)$','Interpreter','Latex','Rotation',90,'FontSize',14)
axis([0 10 0 1])
%fprintf('%2.0f\t %2.0f\t %2.0f\t %2.0f\t %2.0f\t %2.0f\t %2.0f\t %10.8f\t %10.8f\n',G_r,K_r,it,L,N,pr1(N+1),thetar1(N+1));
end
toc
Here is a code for the Blasius equation to give you an idea:
N = 100;
L = 10;
%
% Solve Blasius equation f''' + 0.5*f*f'' = 0, f(0) = f'(0) = 0, f'(1) = 1
% using spectral relaxation method
% u1 = f, u2 = f', u3 = f''
%
maxiter = 40;
tol = 1e-10;
[D,x] = cheb(N);
eta = L/2*(1-x);
u1_old = eta - (1-exp(-eta));
u2_old = 1 - exp(-eta);
u3_old = exp(-eta);
u_old = [u1_old;u2_old;u3_old];
error = 1;
iter = 0;
while error > tol & iter <= maxiter
iter = iter + 1;
A = [-2/L*D,-eye(N+1),zeros(N+1);... % u1' - u2 = 0
zeros(N+1),-2/L*D,-eye(N+1);... % u2' - u3 = 0
zeros(N+1),zeros(N+1),-2/L*D+0.5*diag(u1_old);... % u3' + 0.5*u1*u3 = 0
[1,zeros(1,N)],zeros(1,N+1),zeros(1,N+1);... % u1(0) = 0
zeros(1,N+1),[1,zeros(1,N)],zeros(1,N+1);... % u2(0) = 0
zeros(1,N+1),[zeros(1,N),1],zeros(1,N+1)]; % u2(L) = 1
b = [zeros((N+1),1);zeros(N+1,1);zeros(N+1,1);0;0;1]; % Right-hand sides: 0,0,0,0,0,1
u = A\b;
error = norm(u-u_old,Inf);
u1_old = u(1:N+1);
u2_old = u((N+1)+1:2*(N+1));
u3_old = u(2*(N+1)+1:3*(N+1));
u_old = [u1_old;u2_old;u3_old];
end
iter
error
u1_SRM = u(1:N+1);
u2_SRM = u((N+1)+1:2*(N+1));
u3_SRM = u(2*(N+1)+1:3*(N+1));
%
% Solve Blasius equation f''' + 0.5*f*f'' = 0, f(0) = f'(0) = 0, f'(1) = 1
% using bvp4c
%
xmesh = linspace(0,L,N+1);
solinit = bvpinit(xmesh, [0;0;0]);
odefun = @(x,y)[y(2);y(3);-0.5*y(1)*y(3)];
bcfun = @(ya,yb)[ya(1);ya(2);yb(2)-1];
sol = bvp4c(odefun,bcfun,solinit);
u = deval(sol,eta).';
u1_bvp4c = u(:,1);
u2_bvp4c = u(:,2);
u3_bvp4c = u(:,3);
% Compare results
norm(u1_bvp4c-u1_SRM,Inf)
norm(u2_bvp4c-u2_SRM,Inf)
norm(u3_bvp4c-u3_SRM,Inf)
figure(1)
plot(eta,[u1_SRM,u1_bvp4c])
figure(2)
plot(eta,[u2_SRM,u2_bvp4c])
figure(3)
plot(eta,[u3_SRM,u3_bvp4c])
function [D, x] = cheb(N)
if N == 0
D = 0;
x = 1;
return
end
x = cos(pi * (0:N) / N)';
c = [2; ones(N - 1, 1); 2] .* (-1).^(0:N)';
X = repmat(x, 1, N + 1);
dX = X - X';
D = (c * (1 ./ c)')./(dX + (eye(N + 1)));
D = D - diag(sum(D, 2));
end
Categories
Find more on Chemistry in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
