Error in Finite Difference Method for ODE Code
11 views (last 30 days)
Show older comments
I'm working on an ODE Finite Difference Method. I think I have it down (so I deleted my previous question), but I am getting this error:
Array indices must be positive integers or logical values.
Error in FiniteDifference (line 18)
p(X) = (1 + 0*X); q(X) = (2 + 0*X); r(X) = cos(X);
Here is my code:
function [A,B,W] = FiniteDifference(N)
%equation
%y'' = y' + 2y + cosx, 0<x<pi/2
%y(0) = -0.3, y(pi/2) = -0.1
%breaking it down:
% a = 0; alpha = -0.3; b = pi/2; beta = -0.1;
% p(x) = 1; q(x) = 2; r(x) = cosx;
%to check: plot and residue
A = zeros(1,N+1); B = zeros(1,N+1); C = zeros(1,N+1); D = zeros(1,N+1);
L = zeros(1,N+1); U = zeros(1,N+1); Z = zeros(1,N+1); W = zeros(1,N+1);
aa = 0; alpha = -0.3;
bb = pi/2; beta = -0.1;
h = (bb-aa)/(N+1); %Step Size
X = aa + h;
p(X) = (1 + 0*X); q(X) = (2 + 0*X); r(X) = cos(X);
A(1) = 2+h^2*q(X);
B(1) = -1+0.5*h*p(X);
D(1) = -h^2*r(X)+(1+0.5*h*p(X))*alpha;
for i = 2:N-1
X = aa + i*h;
A(i) = 2+h^2*q(X);
B(i) = -1+0.5*h*p(X);
C(i) = -1-0.5*h*p(X);
D(i) = -h^2*r(X);
end
X = bb-h;
A(N) = 2 + h^2*q(X);
C(N) = -1 - 0.5*h*p(X);
D(N) = -h^2*r(X) + (1 - 0.5*h*p(X))*beta;
L(1) = A(1);
U(1) = B(1)/A(1);
Z(1) = D(1)/L(1);
for i = 2:N-1
L(i) = A(i) - C(i)*U(i-1);
U(i) = B(i)/L(i);
Z(i) = (D(i) - C(i)*Z(i-1))/L(i);
end
L(N) = A(N) - C(N)*U(N-1);
Z(N) = (D(N) - C(N)*Z(N-1))/L(N);
W(N) = Z(N);
for j = 1:N-1
i = N-j;
W(i) = Z(i) - U(i)*W(i+1);
end
i = 0;
for i = 1:N
X = aa + i*h;
end
plot(X,W);
end
0 Comments
Accepted Answer
Torsten
on 30 Apr 2018
Edited: Torsten
on 30 Apr 2018
Maybe you mean
p(X) = @(X)(1 + 0*X); q(X) = @(X)(2 + 0*X); r(X) = @(X)cos(X);
?
Best wishes
Torsten.
4 Comments
Jan
on 2 May 2018
r = @(X)cos(X)? What about r=@cos? Anonymous functions are much slower than function handles. q = @(X)(2 + 0*X)? Looks like: q = 2. Do I oversee something?
More Answers (0)
See Also
Categories
Find more on Mathematics and Optimization in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!