Help me, Error message: "?? Subscript indices must either be real positive integers or logicals"
4 views (last 30 days)
Show older comments
Dear Experts, I'm new in using matlab and I'm not using any programming software before. I wanted to use gaver-stehfest algorithm in inverting a laplace equation.
I made these function codes in three m-files:
Function 1:
function [f,y]=prnfr(s) ome=1; lam=1; f=(ome(1-ome)*s+lam)/((1-ome)*s+lam); y=besselk(0,(sqrt(s*f)))/(s*sqrt(s*f)*besselk(1,(sqrt(s*f))));
Function 2: % % ilt=gavsteh(funname,t,L) % % funname The name of the function to be transformed. % t The transform argument (usually a snapshot of time). % ilt The value of the inverse transform % L number of coefficient ---> depends on computer word length used % (examples: L=8, 10, 12, 14, 16, so on..) % % Wahyu Srigutomo % Physics Department, Bandung Institute of Tech., Indonesia, 2006 % Numerical Inverse Laplace Transform using Gaver-Stehfest method % %Refferences: % 1. Villinger, H., 1985, Solving cylindrical geothermal problems using % Gaver-Stehfest inverse Laplace transform, Geophysics, vol. 50 no. 10 p. % 1581-1587 % 2. Stehfest, H., 1970, Algorithm 368: Numerical inversion of Laplace transform, % Communication of the ACM, vol. 13 no. 1 p. 47-49 % % Simple (and yet rush) examples included in functions fun1 and fun2 with % their comparisons to the exact value (use testgs.m to run the examples) function ilt=gavsteh(funname,t,L)
nn2 = L/2; nn21= nn2+1;
for n = 1:L z = 0.0; for k = floor( ( n + 1 ) / 2 ):min(n,nn2) z = z + ((k^nn2)*factorial(2*k))/ ... (factorial(nn2-k)*factorial(k)*factorial(k-1)* ... factorial(n-k)*factorial(2*k - n)); end v(n)=(-1)^(n+nn2)*z; end
sum = 0.0; ln2_on_t = log(2.0) / t; for n = 1:L p = n * ln2_on_t; sum = sum + v(n) * feval(funname,p); end ilt = sum * ln2_on_t;
Function 3: % Applying Gaver-stehvest to the laplace equation
function y=gs(L); L=input('L = ');
for l=1:L t(l)=l*1; F(l)=gavsteh('prnfr',t(l),L);
end Result1=[t' F'] plot(t,F)
I run gs.m file and I always get this error message: >> gs L = 20 ??? Subscript indices must either be real positive integers or logicals.
Error in ==> prnfr at 4 f=(ome(1-ome)*s+lam)/((1-ome)*s+lam);
Error in ==> gavsteh at 42 sum = sum + v(n) * feval(funname,p);
Error in ==> gs at 9 F(l)=gavsteh('prnfr',t(l),L);
Can anybody help please?
1 Comment
Walter Roberson
on 4 Apr 2011
Please edit the question, select the code, and click on the 'Code {}' button. That will format the code so that it is readable by others.
Answers (2)
Walter Roberson
on 4 Apr 2011
You have ome=1; f=(ome(1-ome)*s+lam)/((1-ome)*s+lam); Did you mean f=(ome*(1-ome)*s+lam)/((1-ome)*s+lam); ?
0 Comments
See Also
Categories
Find more on Transforms 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!