Error in Finite Difference Method for ODE Code

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

 Risposta accettata

Torsten
Torsten il 30 Apr 2018
Modificato: Torsten il 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 Commenti

I completely forgot to add that. Thank you! However, it's still giving me the same error. I'm not sure what else I missed.
Sorry, should have been
p = @(X)(1 + 0*X); q = @(X)(2 + 0*X); r = @(X)cos(X);
Best wishes
Torsten.
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?
Tortsen, thank you so much!!!

Accedi per commentare.

Più risposte (0)

Community Treasure Hunt

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

Start Hunting!

Translated by