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 Kommentar

per isakson
per isakson am 25 Nov. 2020
Write some comments, which explains what the loops are intended to achieve.

Melden Sie sich an, um zu kommentieren.

 Akzeptierte Antwort

VBBV
VBBV am 25 Nov. 2020

0 Stimmen

%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 Kommentare

Naman Singh
Naman Singh am 25 Nov. 2020
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 am 25 Nov. 2020
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 am 25 Nov. 2020
in the same function you are also trying to calculate x vector using the matrix inversion operator \
VBBV
VBBV am 25 Nov. 2020
Also an end statement for the function is required

Melden Sie sich an, um zu kommentieren.

Weitere Antworten (0)

Kategorien

Tags

Gefragt:

am 25 Nov. 2020

Kommentiert:

am 25 Nov. 2020

Community Treasure Hunt

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

Start Hunting!

Translated by