i have error in runge kutta 4 order.

Question

Rashmi 2024-6-13

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2128286-i-have-error-in-runge-kutta-4-order

编辑： Walter Roberson 2024-6-16

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.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

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

Answer 1

Ganesh 2024-6-13

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2128286-i-have-error-in-runge-kutta-4-order#answer_1471491

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 2024-6-16

No S not a function in this line I am forget *.

Torsten 2024-6-16

编辑：Torsten 2024-6-16

在 MATLAB Online 中打开

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));