Clear Filters
Clear Filters

H have a error in FDM.

20 views (last 30 days)
Rashmi
Rashmi on 23 Jun 2024 at 11:21
Answered: Torsten on 23 Jun 2024 at 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.

Accepted Answer

Torsten
Torsten on 23 Jun 2024 at 13:10
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.

More Answers (0)

Categories

Find more on Desktop in Help Center and File Exchange

Tags

Products


Release

R2024a

Community Treasure Hunt

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

Start Hunting!