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 Ansicht (letzte 30 Tage)
Ältere Kommentare anzeigen
% 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 Kommentare
Antworten (1)
Walter Roberson
am 3 Sep. 2022
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 Kommentare
Siehe auch
Kategorien
Mehr zu Function Creation finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!