H have a error in FDM.

Question

Rashmi 2024-6-23，11:21

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2131191-h-have-a-error-in-fdm

回答： Torsten 2024-6-23，13:10

FDMex()
Array indices must be positive integers or logical values.

Error in solution>FDMex/equation (line 61)
       F(i+1) = (F(i) + F(i+2))/2 -H(i)*(h/2)*(F(i+1)-F(i)) +(h^2)*(G(i)^2) - (h^2)*(F(i)^2) + (h^2)*(S/2)*((((i*h)+1)/2)*((F(i+1)-F(i))/h)+F(i)) - (h^2)*(k/2)*(((G(i+1)-G(i))/h)^2 - ((F(i+1)-F(i))/h)^2+2*F((F(i) -2*F(i+1) + F(i+2))/h^2)); 

Error in solution>FDMex (line 34)
    [H1,F1,G1,theta1,phi1] = equation(H,F,G,theta,phi,N,h); 
function FDMex() 
N = 100; 
lgth = 1.0; 
h = lgth/N; 
eta = 0:h:lgth; 
Ac = 0.0001; 
S = 0.2; 
k = 0.1; 
Pr = 1.0; 
Sc = 1.2; 
alpha1 = 0.4; 
alpha2 = 0; 
zeta = 0.3; 
gamma = 0.3; 
omega = 0.4; 
fw = 0.2;
F = zeros(N+2, 1); 
G = zeros(N+2, 1); 
theta = zeros(N+2, 1); 
phi = zeros(N+2, 1); 
H = zeros(N+2, 1);
F(1) = 0; 
G(1) = omega; 
theta(1) = 1; 
phi(1) = 1; 
H(1) = S + fw / Sc * (phi(2) - phi(1)) / h^2; 
F(N+2) = 0; 
G(N+2) = 0; 
theta(N+2) = 0; 
phi(N+2) = 0;
c = 1.0; 
while(c>0) 
    [H1,F1,G1,theta1,phi1] = equation(H,F,G,theta,phi,N,h); 
    c = 0.0; 
    for i = 2:N-1 
        if (abs((H(i)-H1(i)), (F(i) - F1(i)), (G(i) - G1(i)),(theta(i) - theta1(i)), (phi(i) - phi1(i)))>Ac) 
            c = c+1; 
            break 
        end 
    end 
    H = H1; 
    F = F1; 
    G = G1; 
    theta = theta1; 
    phi = phi1; 
end 
disp('Hence solutions = :' ); 
H2(1 : N+2) = H; 
F2(1 : N+2) = F; 
G2(1 : N+2) = G; 
theta2(1 : N+2) = theta; 
phi2(1 : N+2) = phi; 
eta = 0:h:lgth;
figure(1) 
plot(eta,H2,'*r') 
hold on
function [H1,F1,G1,theta1,phi1] = equation(H,F,G,theta,phi,N,h) 
   for i = 2:N-1 
       H(i+1) = H(i) - h*2*F(i); 
       F(i+1) = (F(i) + F(i+2))/2 -H(i)*(h/2)*(F(i+1)-F(i)) +(h^2)*(G(i)^2) - (h^2)*(F(i)^2) + (h^2)*(S/2)*((((i*h)+1)/2)*((F(i+1)-F(i))/h)+F(i)) - (h^2)*(k/2)*(((G(i+1)-G(i))/h)^2 - ((F(i+1)-F(i))/h)^2+2*F((F(i) -2*F(i+1) + F(i+2))/h^2)); 
       G(i+1) = 2*G(i) - G(i+2) -(h/2)*H(i)*(G(i+1)-G(i+2)) + 2*F(i)*G(i) - (h^2)*S*(((i*h+1)/2)*((G(i+1)-G(i+2))/2*h)+G(i)) + (h^2)*2*k*(F(i)*((G(i+1) -2*G(i) + G(i+2))/h^2) - ((F(i+1)-F(i+2))/2*h)*((G(i+1)-G(i+2))/2*h)); 
       theta(i+1) = 2*theta(i) - theta(i+2) + Pr*(h/2)*(theta(i+1)-theta(i+2)) - Pr*(h^2)*S*(((i*h+1)/2)*((theta(i+1)-theta(i+2))/2*h)+alpha1*theta(i)) + zeta*((theta(i+1) -2*theta(i) + theta(i+2))/h^2 - 2*F(i)*G(i)*((theta(i+1)-theta(i+2))/2*h)); phi(i+1) = 2*phi(i) - phi(i+2) + Sc*(h/2)*(theta(i+1)-theta(i+2)) - Sc*(h^2)*S*(((i*h+1)/2)*((phi(i+1)-theta(i+2))/2*h)+alpha2*phi(i)) + Sc*(h^2)*gamma*phi(i); 
   end 
   H1(1) = H(1); 
   F1(1) = F(1); 
   G1(1) = G(1); 
   theta1(1) = theta(1); 
   phi1(1) = phi(1);
   F1(N+2) = F(N+2); 
   G1(N+2) = G(N+2); 
   theta1(N+2) = theta(N+2); 
   phi1(N+2) = phi(N+2); 
end 
end 

AArray indices must be positive integers or logical values.

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

请先登录，再进行评论。

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

Answer 1

Torsten 2024-6-23，13:10

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2131191-h-have-a-error-in-fdm#answer_1475911

在 MATLAB Online 中打开

       F(i+1) = (F(i) + F(i+2))/2 -H(i)*(h/2)*(F(i+1)-F(i)) +(h^2)*(G(i)^2) - (h^2)*(F(i)^2) + (h^2)*(S/2)*((((i*h)+1)/2)*((F(i+1)-F(i))/h)+F(i)) - (h^2)*(k/2)*(((G(i+1)-G(i))/h)^2 - ((F(i+1)-F(i))/h)^2+2*F((F(i) -2*F(i+1) + F(i+2))/h^2)); 
       G(i+1) = 2*G(i) - G(i+2) -(h/2)*H(i)*(G(i+1)-G(i+2)) + 2*F(i)*G(i) - (h^2)*S*(((i*h+1)/2)*((G(i+1)-G(i+2))/2*h)+G(i)) + (h^2)*2*k*(F(i)*((G(i+1) -2*G(i) + G(i+2))/h^2) - ((F(i+1)-F(i+2))/2*h)*((G(i+1)-G(i+2))/2*h)); 
       theta(i+1) = 2*theta(i) - theta(i+2) + Pr*(h/2)*(theta(i+1)-theta(i+2)) - Pr*(h^2)*S*(((i*h+1)/2)*((theta(i+1)-theta(i+2))/2*h)+alpha1*theta(i)) + zeta*((theta(i+1) -2*theta(i) + theta(i+2))/h^2 - 2*F(i)*G(i)*((theta(i+1)-theta(i+2))/2*h)); phi(i+1) = 2*phi(i) - phi(i+2) + Sc*(h/2)*(theta(i+1)-theta(i+2)) - Sc*(h^2)*S*(((i*h+1)/2)*((phi(i+1)-theta(i+2))/2*h)+alpha2*phi(i)) + Sc*(h^2)*gamma*phi(i); 
       

Line 1:

You forgot the loop index i you refer to:

2*F((F(i) -2*F(i+1) + F(i+2))/h^2));

Line 1,2 and 3:

You reference array elements of F, G and theta on the right-hand sides that are not yet defined:

If you define F(i+1), G(i+1) and theta(i+1), the expressions on the right-hand side can only reference F(k), G(k) and theta(k) with k < i+1.

Lines 2 and 3:

If you use /2*h, you want to divide by 2*h, but you divide by 2 and multiply the resulting expression by h. Use /(2*h) instead.