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!