i have error in runge kutta 4 order.
3 visualizzazioni (ultimi 30 giorni)
Mostra commenti meno recenti
Alpha1=1;
Alpha2=2;
S=1;
omega=1
gamma=1;
Sc=1;
Pr=1.0;
f_omega=0.4;
zeta=0;
k = 0.1;
h=0.05; % step size
eta = 0:h:100; % Calculates upto y(3)
H = zeros(1,length(eta));
g = zeros(1,length(eta));
F = zeros(1,length(eta));
z = zeros(1,length(eta));
G = zeros(1,length(eta));
x = zeros(1,length(eta));
theta = zeros(1,length(eta));
y = zeros(1,length(eta));
phi = zeros(1,length(eta));
F(1) = 0;
G(1) = omega;
theta(1) = 1;
phi(1) = 1;
H(1) = S +f_omega / Sc * y(1) ;
u1 = @(eta, H) (-2*F);
u2 = @(eta,F) g;
u3 = @(eta,F,g) (g*H - G.^2 +F.^2 - S(((eta+1)/2)*g + F) + k(z.^2 -g.^2)) / (1-2*F*k);
u4 = @(eta,G) z;
u5 = @(eta,G,z) ((H*z + 2*F*G - S(((eta+1)/2)*z + G) - 2*k*g*z)/(1-2*F*k));
u6 = @(eta, theta) x;
u7 = @(eta, theta , x) (Pr*(H*x + S(((rta+1)/2)*x + alpha1*theta) - 2*zeta*F*H*x)/(1-Pr*zeta));
u8 = @(eta, phi) y ;
u9 = @(eta, phi, y) (Sc*(H*y - S(((eta+1)/2)*y + alpha2*phi) +gamma*phi));
for i=1:(length(eta)-1) % calculation loop
k_1 =u1(eta(i),H(i));
l_1 =u2(eta(i),F(i));
m_1 =u3(eta(i),F(i),g(i));
n_1 =u4(eta(i),G(i));
o_1 =u5(eta(i),G(i),z(i));
p_1 =u6(eta(i),theta(i));
q_1 =u7(eta(i),theta(i),x(i));
r_1 =u8(eta(i),phi(i));
s_1 =u9(eta(i),phi(i),y(i));
k_2 = u1(eta(i)+0.5*h,H(i)+0.5*h*k_1);
l_2 = u2(eta(i)+0.5*h,F(i)+0.5*h*l_1);
m_2 = u3(eta(i)+0.5*h,F(i)+0.5*h*l_1,g(i)+0.5*h*m_1);
n_2 = u4(eta(i)+0.5*h,G(i)+0.5*h*n_1);
o_2 = u5(eta(i)+0.5*h,G(i)+0.5*h*n_1,z(i)+0.5*h*o_1);
p_2 = u6(eta(i)+0.5*h,theta(i)+0.5*h*p_1);
q_2 = u7(eta(i)+0.5*h,theta(i)+0.5*h*p_1,x(i)+0.5*h*q_1);
r_2 = u8(eta(i)+0.5*h,phi(i)+0.5*h*r_1);
s_2 = u9(eta(i)+0.5*h,phi(i)+0.5*h*r_1,x(i)+0.5*h*s_1);
k_3 = u1(eta(i)+0.5*h,H(i)+0.5*h*k_2);
l_3 = u2(eta(i)+0.5*h,F(i)+0.5*h*l_2);
m_3 = u3(eta(i)+0.5*h,F(i)+0.5*h*l_2,g(i)+0.5*h*m_2);
n_3 = u4(eta(i)+0.5*h,G(i)+0.5*h*n_2);
o_3 = u5(eta(i)+0.5*h,G(i)+0.5*h*n_2,z(i)+0.5*h*o_2);
p_3 = u6(eta(i)+0.5*h,theta(i)+0.5*h*p_2);
q_3 = u7(eta(i)+0.5*h,theta(i)+0.5*h*p_2,x(i)+0.5*h*q_2);
r_3 = u8(eta(i)+0.5*h,phi(i)+0.5*h*r_2);
s_3 = u9(eta(i)+0.5*h,phi(i)+0.5*h*r_2,x(i)+0.5*h*s_2);
k_4 = u1(eta(i)+h,H(i)+h*k_3);
l_4 = u2(eta(i)+h,F(i)+h*l_3);
m_4 = u3(eta(i)+h,F(i)+h*l_3,g(i)+h*m_3);
n_4 = u4(eta(i)+h,G(i)+h*n_3);
o_4 = u5(eta(i)+h,G(i)+h*n_3,z(i)+h*o_3);
p_4 = u6(eta(i)+h,theta(i)+0.5*h*p_2);
q_4 = u7(eta(i)+h,theta(i)+h*p_3,x(i)+h*q_3);
r_4 = u8(eta(i)+h,phi(i)+h*r_3);
s_4 = u9(eta(i)+h,phi(i)+h*r_3,x(i)+h*s_3);
H(i+1) = H(i) + (1/6)*(k_1+2*k_2+2*k_3+k_4)*h;
g(i+1) = g(i) + (1/6)*(l_1+2*l_2+2*l_3+l_4)*h;
F(i+1) = F(i) + (1/6)*(m_1+2*m_2+2*m_3+m_4)*h;
z(i+1) = z(i) + (1/6)*(n_1+2*n_2+2*n_3+n_4)*h;
G(i+1) = G(i) + (1/6)*(o_1+2*o_2+2*o_3+o_4)*h;
x(i+1) = x(i) + (1/6)*(p_1+2*p_2+2*p_3+p_4)*h;
theta(i+1) = theta(i) + (1/6)*(q_1+2*q_2+2*q_3+q_4)*h;
y(i+1) = y(i) + (1/6)*(r_1+2*r_2+2*r_3+r_4)*h;
phi(i+1) = phi(i) + (1/6)*(s_1+2*s_2+2*s_3+s_4)*h;
end
figure();
subplot(3,2,1);
plot(eta, F);
title('F vs \eta');
xlabel('\eta');
ylabel('F');
subplot(3,2,2);
plot(eta, G);
title('G vs \eta');
xlabel('\eta');
ylabel('G');
subplot(3,2,3);
plot(eta, theta);
title('\theta vs \eta');
xlabel('\eta');
ylabel('\theta');
subplot(3,2,4);
plot(eta, phi);
title('\phi vs \eta');
xlabel('\eta');
ylabel('\phi');
subplot(3,2,5);
plot(eta, H);
title('H vs \eta');
xlabel('\eta');
ylabel('H');
error
Array indices must be positive integers or logical values.
0 Commenti
Risposte (1)
Ganesh
il 13 Giu 2024
The error occurs as you are trying to index the variable "S" while calling function "u3" whereas you have declared S to be an integer value of 1.
u3 = @(eta,F,g) (g*H - G.^2 +F.^2 - S(((eta+1)/2)*g + F) + k(z.^2 -g.^2)) / (1-2*F*k);
2 Commenti
Torsten
il 16 Giu 2024
Modificato: Torsten
il 16 Giu 2024
u1 = @(eta, H) (-2*F);
u2 = @(eta,F) g;
u3 = @(eta,F,g) (g*H - G.^2 +F.^2 - S(((eta+1)/2)*g + F) + k(z.^2 -g.^2)) / (1-2*F*k);
u4 = @(eta,G) z;
u5 = @(eta,G,z) ((H*z + 2*F*G - S(((eta+1)/2)*z + G) - 2*k*g*z)/(1-2*F*k));
u6 = @(eta, theta) x;
u7 = @(eta, theta , x) (Pr*(H*x + S(((rta+1)/2)*x + alpha1*theta) - 2*zeta*F*H*x)/(1-Pr*zeta));
u8 = @(eta, phi) y ;
u9 = @(eta, phi, y) (Sc*(H*y - S(((eta+1)/2)*y + alpha2*phi) +gamma*phi));
Why don't you define u1 to depend on F ? You only make it depend on eta and H.
Why don't you define u2 to depend on g ?
Why don't you define u3 to depend on G and H ?
...
Vedere anche
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!