H have a error in FDM.

조회 수: 4 (최근 30일)
Rashmi
Rashmi 2024년 6월 23일
답변: Torsten 2024년 6월 23일
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.

채택된 답변

Torsten
Torsten 2024년 6월 23일
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.

추가 답변 (0개)

카테고리

Help CenterFile 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