function [x, y1, y2, y3, y4 ] = FDEx4_1(N,kn,F)
clc
mm = 4;
N = 200;
h = 1/N;% mesh size
kn = 0.5;%Reynolds number
F = 0.23;
x = (0:h:1)'; % nodes
M = zeros(mm*(N+1),mm*(N+1));
b = zeros(mm*(N+1),1);
Done = 0;
Max_iter = 10;
y0 = zeros(mm*(N+1),1);
for i = 1:N+1
y0(mm*(i-1)+1) = 0;
y0(mm*(i-1)+2) = 0;
y0(mm*(i-1)+3) = 0;
y0(mm*(i-1)+4) = 0;
end
y1 = zeros(N+1,1);
y2 = zeros(N+1,1);
y3 = zeros(N+1,1);
y4 = zeros(N+1,1);
for i = 1:N+1
y1(i) = x(mm*(i-1)+1);
y2(i) = x(mm*(i-1)+2);
y3(i) = x(mm*(i-1)+3);
y4(i) = x(mm*(i-1)+4);
end
iter = 0;
while (~Done)
for i = 1:N
xip1 = 0.5*(x(i)+x(i+1));
yip1 = 0.5*(y0(mm*(i-1)+1:mm*i)+y0(mm*i+1:mm*(i+1)));
Aip1= [0,-1./kn, 0, 0;
0, 0, -F./(yip1(3)).^2, 0;
0, 0, 0, -4./(15*kn);
0, (2*yip1(2))./(kn), 0, 0];
Qip1 = [-(yip1(2))./kn; F./yip1(3); (4*yip1(4))./(15*kn); ((yip1(2)).^2)./kn];
fip1 = Qip1-Aip1*yip1;
M(mm*(i-1)+1:mm*i,mm*(i-1)+1:mm*i) = -(0.5*Aip1+eye(mm)/h);
M(mm*(i-1)+1:mm*i,mm*i+1:mm*(i+1)) = (-0.5*Aip1+eye(mm)/h);
b(mm*(i-1)+1:mm*i) = fip1;
end
Ba = [sqrt(2./pi*y3(0)), 1, -(y1(0))*(1./sqrt(pi*2)*(y3(0).^1.5)), 0;
0, 0, 0, 0;
0, 0, (sqrt(2)*(1./(y3(0)).^1.5)*sqrt(pi))*(y3(0) + 1), 1;
0, 0, 0, 0];
Bb = [0, 0, 0, 0;
-sqrt(2./(pi*y3(1))), 1, y1(1)*(1./(y3(0)).^1.5*sqrt(2*pi)), 0;
0, 0, 0, 0;
0, 0, -(sqrt(2)*(1./(y3(0)).^1.5)*sqrt(pi))*(y3(1) + 1), 1];
M(mm*((N+1)-1)+1:mm*(N+1),mm*((N+1)-1)+1:mm*(N+1)) = Bb;
M(mm*((N+1)-1)+1:mm*(N+1),1:mm) = Ba;
alpha = [0; 0; 0; 0];
b(mm*((N+1)-1)+1:mm*(N+1)) = alpha;
x = M\b;
error = max(abs(x-y0));% max |y_i-y_i^*|
iter = iter+1;
disp(['Error in step ',num2str(iter),' is : ',num2str(error)])
y0 = x;
Done = (iter>=Max_iter)||(error<0.000001);
end
It says
Index exceeds the number of array elements (201).
Error in FDEx4_1 (line 31)
y2(i) = x(mm*(i-1)+2);
How do i fix this?

댓글 수: 1

per isakson
per isakson 2020년 11월 25일
Write some comments, which explains what the loops are intended to achieve.

댓글을 달려면 로그인하십시오.

 채택된 답변

VBBV
VBBV 2020년 11월 25일

0 개 추천

%if true
x = 0:h:1;% size has 201 elements
Define this according to function calculations e.g.
x = 0:1/((N+1)*mm):1;% size of 805

댓글 수: 4

Naman Singh
Naman Singh 2020년 11월 25일
Ba = [sqrt(2./pi*y3(0)), 1, -(y1(0))*(1./sqrt(pi*2)*(y3(0).^1.5)), 0;
In this line it gives an error of
Array indices must be positive integers or logical values.
VBBV
VBBV 2020년 11월 25일
Matlab does not use zero or decimal indexing for matrix definitions. Since y1 and y3 are matrices using 0 as index which is undefined in matlab, it produces error.
You can set it to y1(1) or y3(1) if you want to use only first element of those matrix
VBBV
VBBV 2020년 11월 25일
in the same function you are also trying to calculate x vector using the matrix inversion operator \
VBBV
VBBV 2020년 11월 25일
Also an end statement for the function is required

댓글을 달려면 로그인하십시오.

추가 답변 (0개)

카테고리

도움말 센터File Exchange에서 Matrix Indexing에 대해 자세히 알아보기

태그

질문:

2020년 11월 25일

댓글:

2020년 11월 25일

Community Treasure Hunt

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

Start Hunting!

Translated by