Error in Finite Difference Method for ODE Code

13 次查看（过去 30 天）
Mel Hernandez
Mel Hernandez 2018-4-30
评论： Mel Hernandez 2018-5-2
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

请先登录，再回答此问题。

采纳的回答

Torsten
Torsten 2018-4-30
编辑：Torsten 2018-4-30
Maybe you mean
p(X) = @(X)(1 + 0*X); q(X) = @(X)(2 + 0*X); r(X) = @(X)cos(X);
?
Best wishes
Torsten.
  4 个评论
显示 2更早的评论
Jan
Jan 2018-5-2
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?

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Mathematics and Optimization 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by