フィルターのクリア

H have a error in FDM.

6 ビュー (過去 30 日間)
Rashmi
Rashmi 2024 年 6 月 23 日 11:21
回答済み: 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.

採用された回答

Torsten
Torsten 2024 年 6 月 23 日 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.

その他の回答 (0 件)

カテゴリ

Help Center および File ExchangeDesktop についてさらに検索

タグ

製品


リリース

R2024a

Community Treasure Hunt

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

Start Hunting!

Translated by