Unable to perform assignment because the size of the left side is 1-by-3 and the size of the right side is 1-by-2.
1 次查看(过去 30 天)
显示 更早的评论
% system of first order differential equation
% y1' = y2;
% y2' = y3;
% y3' = (1-phi).^2.5*[(1-phi)+phi*(rows/rowp)*(y1*y3-y1.^2)-M*y2+lamda*y4*((1-phi)+phi*(rowbetas/rowbetaf))]
% y4' = y5;
% y5' = -(3*N*pr*kf*[(1-phi)+phi*(rowcps/rowcpf)]/(3*N+4)*knf)*y1*y5;
%with boundary condition y1(0) = 0, y2(0)=1, y3(0)=0.078, y4(0) = 1, y5(0)=?
phi = 0.01;
lamda = 1;
M = 2;
Ec = 2;
N = 1.5;
pr = 6.2;
%copper
rows =8933 ;
ks = 401;
cps =385 ;
rowcps=8933*385;
rowbetas = 8933*1.67*10.^-5;
betas = 1.67*10.^-5;
%alumina
% rows =3970 ;
% ks =40 ;
% cps =765;
%betas = 0.85*10.^-5;
%water
rowf = 997.1;
kf = 0.613;
cpf = 4179;
betaf = 21*10.^-5;
rowbetaf = 997.1*21*10.^-5;
rowcpf = 997.1*4179;
% constants
A1 = (1-phi)+phi*(rows/rowf);
A2 = (1-phi)+phi*(rowbetas/rowbetaf);
A3 = (1-phi)+phi*(rowcps/rowcpf);
A = (1-phi).^2.5;
A4 = -(3*N*pr*kf/(3*N+4)*knf);
knf = kf*(ks+2*kf-2*phi*(kf-ks)/ks+2*kf+phi*(kf-ks));
%% solve ODE-IVP using RK-4 Method
tic;
% range of independent variable η
eta_0 = 0;
eta_max = 5;
% Initial conditions
f0 = 0;
g0 = 1;
h0 = -0.74878322
l0 = 1;
m0 = 0.2356
% Step value Δη
d_eta = 0.001;
N = (eta_max - eta_0)/ d_eta;
err= 1;
%% Initializing solution
eta = eta_0 : d_eta : eta_max; % X coordinate for plot function
while err >10^-9
F = zeros(N+1,3); % f = F(N+1,1) ; g = F(N+1,2) ; h = F(N+1,3) ;
F(1,:) = [f0 g0 h0];
G = zeros(N+1,2);
G(1,:) = [l0 m0]; %l = G(N+1,1); m = G(N+1,2); j let us choose artbitrary for two values
%% Trail using RK-4 Method
for i = 1:N
k11(i,:) = d_eta*[G(i,2) A4*A3*F(i,1)*G(i,2)];
k1(i,:) = d_eta* [F(i,2) F(i,3) A*(A1*(F(i,1)*F(i,3)-F(i,2).^2)-M*F(i,2)+lamda*A2*G(i,1))] ;
k12(i,:) = d_eta*[G(i,2)+k11(i,1)/2 A4*A3*F(i,1)+k1(i,1)/2*G(i,2)+k11(i,2)/2];
k2(i,:) = d_eta* [F(i,2)+k1(i,2)/2 F(i,3)+k1(i,3)/2 A*(A1*(F(i,1)+k1(i,1)/2*F(i,3)+k1(i,3)/2-F(i,2).^2+k1(i,2)/2)-M*F(i,2)+k1(i,2)/2+lamda*A2*G(i,1)+k11(i,1)) ];
k13(i,:) = d_eta*[G(i,2)+k12(i,1)/2 A4*A3*F(i,1)+k2(i,1)/2*G(i,2)+k12(i,2)/2];
k3(i,:) = d_eta* [F(i,2)+k2(i,2)/2 F(i,3)+k2(i,3)/2 A*(A1*(F(i,1)+k2(i,1)/2*F(i,3)+k2(i,3)/2-F(i,2).^2+k2(i,2)/2)-M*F(i,2)+k2(i,2)/2+lamda*A2*G(i,1)+k12(i,1)/2) ];
k14(i,:) = d_eta*[G(i,2)+k13(i,1) A4*A3*F(i,1)+k3(i,1)*G(i,2)+k13(i,2)];
k4(i,:) = d_eta* [F(i,2)+k3(i,2) F(i,3)+k3(i,3) A*(A1*(F(i,1)+k3(i,1)*F(i,3)+k3(i,3)-F(i,2).^2+k3(i,2))-M*F(i,2)+k3(i,2)+lamda*A2*G(i,1)+k13(i,1))];
G(i+1,:) = G(i,:)+(1/6)*(k11(i,:)+2*k12(i,:)+2*k13(i,:)+k14(i,:));
F(i+1,:) = F(i,:)+(1/6)*(k1(i,:)+2*k2(i,:)+2*k3(i,:)+k4(i,:));
end
end
%% plotting
f1=figure;
f1.Units = 'normalized';
f1.Position = [0.1 0.1 0.8 0.6];
%plot(eta',F,'LineWidth',1);
plot(eta',G,'Linewidth',1);
0 个评论
回答(1 个)
Walter Roberson
2022-9-3
d_eta = 0.001;
N = (eta_max - eta_0)/ d_eta;
You use N as a size, so you are assuming that it is an integer. However, 0.001 is not exactly representatable in IEEE 754 double precision floating point, so dividing by 0.001 is not certain to give you back an integer, and is not exactly the same as multiplying by 1000.
0 个评论
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Function Creation 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!