I have error in FDM 2

Question

Rashmi 2024-6-23，13:49

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2131206-i-have-error-in-fdm-2

评论： Walter Roberson 2024-6-24，2:37

function FDMex() N = 100;

h = 0.1; eta = 0:h:N; 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 j = 1:N+1 if (abs(H(j)-H1(j))>Ac); (abs(F(j)-F1(j))>Ac);(abs(G(j)-G1(j))>Ac);(abs(theta(1)-theta1(1))>Ac);(abs(phi(j)-phi1(j))>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:0.1:N;

figure(1) plot(eta,H2,'*r') hold on

function [H1, F1, G1, theta1, phi1] = equation(H,F,G,theta,phi,N,h) for i = 1: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(i)*((F(i) -2*F(i+1) + F(i+2))/h^2)); G(i+1) = (G(i) + G(i+2))/2 -(h/2)*H(i)*(G(i+1)-G(i)) - (h^2)*F(i)*G(i) + (h^2)*S*(((i*h+1)/2)*((G(i+1)-G(i))/h)+G(i)) - (h^2)*k*(F(i)*((G(i) -2*G(i+1) + G(i+2))/h^2) - ((F(i+1)-F(i))/h)*((G(i+1)-G(i))/h)); theta(i+1) = (theta(i) + theta(i+2))/2 - Pr*(h/2)*(theta(i+1)-theta(i))*H(i) + Pr*(h^2)*(S/2)*(((i*h+1)/2)*((theta(i+1)-theta(i))/h)+alpha1*theta(i)) - zeta*Pr*((h^2)/2)*((theta(i) -2*theta(i+1) + theta(i+2))/h^2 - 2*F(i)*H(i)*((theta(i+1)-theta(i))/h)); phi(i+1) = (phi(i) + phi(i+2))/2 + Sc*(h/2)*H(i)*(theta(i+1)-theta(i)) - Sc*(h^2)*(S/2)*(((i*h+1)/2)*((phi(i+1)-phi(i))/2*h)+alpha2*phi(i)) + (Sc/2)*(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 Index exceeds the number of array elements. Index must not exceed 1. Line 39

6 个评论
显示 4更早的评论隐藏 4更早的评论

Rashmi 2024-6-23，18:08

function FDMex() N = 100;

h = 0.1; eta = 0:h:N; 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 j = 1:N+1 if (abs(H(j)-H1(j))>Ac); (abs(F(j)-F1(j))>Ac);(abs(G(j)-G1(j))>Ac);(abs(theta(1)-theta1(1))>Ac);(abs(phi(j)-phi1(j))>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:0.1:N;

figure(1) plot(eta,H2,'*r') hold on

function [H1, F1, G1, theta1, phi1] = equation(H,F,G,theta,phi,N,h) for i = 1: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(i)*((F(i) -2*F(i+1) + F(i+2))/h^2)); G(i+1) = (G(i) + G(i+2))/2 -(h/2)*H(i)*(G(i+1)-G(i)) - (h^2)*F(i)*G(i) + (h^2)*S*(((i*h+1)/2)*((G(i+1)-G(i))/h)+G(i)) - (h^2)*k*(F(i)*((G(i) -2*G(i+1) + G(i+2))/h^2) - ((F(i+1)-F(i))/h)*((G(i+1)-G(i))/h)); theta(i+1) = (theta(i) + theta(i+2))/2 - Pr*(h/2)*(theta(i+1)-theta(i))*H(i) + Pr*(h^2)*(S/2)*(((i*h+1)/2)*((theta(i+1)-theta(i))/h)+alpha1*theta(i)) - zeta*Pr*((h^2)/2)*((theta(i) -2*theta(i+1) + theta(i+2))/h^2 - 2*F(i)*H(i)*((theta(i+1)-theta(i))/h)); phi(i+1) = (phi(i) + phi(i+2))/2 + Sc*(h/2)*H(i)*(theta(i+1)-theta(i)) - Sc*(h^2)*(S/2)*(((i*h+1)/2)*((phi(i+1)-phi(i))/2*h)+alpha2*phi(i)) + (Sc/2)*(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

Torsten 2024-6-23，18:13

编辑：Torsten 2024-6-23，18:25

Even without formatting, I see that the main problem is not solved, namely that you want to compute values F(i+1), G(i+1) and theta(i+1) with values F(k), G(k) and theta(k) with k >= i+1. That doesn't make sense.

And the division /(2*h) instead of /2*h has also not been corrected.

You should include a mathematical description of your problem - maybe someone can give you a hint for a valid solution.

Walter Roberson 2024-6-24，2:37

@Rashmi

First press the button marked '>' in the CODE section of the edit toolstrip. This will open a code edit region. Then copy and paste the code from your editor into the code edit region that was opened.

请先登录，再进行评论。

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