function [x, y1, y2, y3, y4 ] = FDEx4_1(N,kn,F)
clc
mm = 4;
N = 200;
h = 1/N;
kn = 0.5;
F = 0.23;
x = (0:h:1)';
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));
iter = iter+1;
disp(['Error in step ',num2str(iter),' is : ',num2str(error)])
y0 = x;
Done = (iter>=Max_iter)||(error<0.000001);
end
1 Comment
Direct link to this comment
https://se.mathworks.com/matlabcentral/answers/661708-index-exceeds-the-number-of-array-elements#comment_1159458
Direct link to this comment
https://se.mathworks.com/matlabcentral/answers/661708-index-exceeds-the-number-of-array-elements#comment_1159458
Sign in to comment.