i have error in runge kutta 4 order.

2 次查看（过去 30 天）
Rashmi
Rashmi 2024-6-13
编辑： Walter Roberson 2024-6-16
Ran in:
Alpha1=1;
Alpha2=2;
S=1;
omega=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
Array indices must be positive integers or logical values.

Error in solution>@(eta,F,g)(g*H-G.^2+F.^2-S(((eta+1)/2)*g+F)+k(z.^2-g.^2))/(1-2*F*k) (line 31)
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);
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.

请先登录，再回答此问题。

回答（1 个）

Ganesh
Ganesh 2024-6-13
Hi @Rashmi,
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 个评论
Rashmi
Rashmi 2024-6-16
No S not a function in this line I am forget *.
Torsten
Torsten 2024-6-16
编辑：Torsten 2024-6-16
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 ?
...

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Loops and Conditional Statements 的更多信息

标签

产品


版本

R2024a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by