polyval error with recurrence relation

T
T
2013 10 6
2 个回答
回答已采纳
更新时间：2013 10 11
9 次查看（30 天）

0 个投票

function c = relation( n )
if n==0
c = 1;
elseif n==1
syms x
c = [x 0];
else
syms x
c = ((2*n-1)*x*[relation(n-1),0] - (n-1)*[0,0,relation(n-2)])/n;
sum(c)
end
n=8
m=n+1;
xval=linspace(-1,1,5000)
I=trapsum(polyval(relation(n),xval).*polyval(relation(m),xval),-1,1);
Gives the following output:
??? Error using ==> polyval
Inputs to polyval must be floats, namely single or double.
Error in ==> polyval at 63
y = zeros(siz_x, superiorfloat(x,p));
What needs to be done to ensure that the output for the recurrence relation is a double?
I mean, in the function, we have syms x but this is to acquire the polynomial in analytic form. It's not necessary in this case but I need to plot it later on. The main issue I have is evaluating the integral using trapsum but I would need to do both; plot for various n, evaluate the integral.

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

 采纳的回答

Walter Roberson
Walter Roberson 2013-10-6

0 个投票

In your code for relation(), why do you have
sum(c)
which is just displaying the sum rather than returning the sum ?
After you build the relation symbolically, use
RelationM = relation(m);
RelationN = relation(n);
funM = matlabFunction(RelationM);
funN = matlabFunction(RelationN);
Y = funM(xval) .* funN(xval);
trapsum(Y, -1, 1)
However, I do not know what your trapsum routine does: it is not a routine in any Mathworks toolbox, and the obvious trapsum in the File Exchange uses a different calling sequence.

14 个评论
显示 12更早的评论 隐藏 12更早的评论

T
T 2013-10-6
Oh. What is the command to return the sum?
I was having difficulties returning the generated polynomial so i just output the result and plot the polynomial manually.
T
T 2013-10-6
I never heard of matlabFunction...interesting...
function c = relation( n )
if n==0
c = 1;
elseif n==1
syms x
c = [x 0];
else
syms x
c = ((2*n-1)*x*[relation(n-1),0] - (n-1)*[0,0,relation(n-2)])/n;
sum(c)
end
n=8
m=n+1;
xval=linspace(-1,1,5000)
RelationM = relation(m);
RelationN = relation(n);
funM = matlabFunction(RelationM);
funN = matlabFunction(RelationN);
Y = funM(xval) .* funN(xval);
trapsum(Y, -1, 1)
Gives me the following:
??? Error using ==> times
Matrix dimensions must agree.
T
T 2013-10-6
Also suppose one wanted to plot:
n=8
m=n+1;
x =[-5:0.05:5];
y = funM
figure(1)
plot(x,y)
Then I receive the message as well:
??? Error using ==> plot
Conversion to double from function_handle is not possible.
Walter Roberson
Walter Roberson 2013-10-6
编辑：Walter Roberson 2013-10-6
To return the sum,
c = sum(c);
In your small test,
plot(x, funM(x))
T
T 2013-10-7
编辑：T 2013-10-7
Adding c = sum(c); Gives the following:
function c = relation( n )
if n==0
c = 1;
elseif n==1
syms x
c = [x 0];
else
syms x
c = ((2*n-1)*x*[relation(n-1),0] - (n-1)*[0,0,relation(n-2)])/n;
c = sum(c)
end
n=3
m=n+1;
x =[-5:0.05:5];
RelationM = relation(m);
RelationN = relation(n);
funM = matlabFunction(RelationM);
funN = matlabFunction(RelationN);
plot(x, funM(x))
??? Error using ==> mupadmex
Error in MuPAD command: Array sizes must match.
Error in ==> sym.minus at 14
X = mupadmex('mllib::zip',A.s,B.s,'_subtract');
Error in ==> relation at 9
c = ((2*n-1)*x*[relation(n-1),0] - (n-1)*[0,0,relation(n-2)])/n;
Error in ==> relation at 9
c = ((2*n-1)*x*[relation(n-1),0] - (n-1)*[0,0,relation(n-2)])/n;
Error in ==> Q922 at 4
RelationM = relation(m);
Then I get the same message if I replace c = sum(c) with just sum(c) at:
??? Error using ==> times
Matrix dimensions must agree.
Error in ==> Q922 at 8
Y = funM(xval) .* funN(xval);
Walter Roberson
Walter Roberson 2013-10-7
Break your expression for c into two pieces, one for relation(n-1) and the other for relation(n-2). Check that their lengths are equal before trying to subtract them. If they are not equal, dump out a warning message. Put a breakpoint in where the warning message would be given, and run the program, and when it stops, examine the variables.
T
T 2013-10-7
Right:
size(funM(xval)) gives 50003
size(funN(xval)) gives 50002
The point was to integrate P_i * P_i+1 on [-1,1]
Walter Roberson
Walter Roberson 2013-10-7
Yes, the two will be different by one element, but the code around the first call adds a 0 after the first vector, and the code around the second call adds two leading 0's before the second vector, so the two should end up the same length.
Go for the full surrounding (2*n-1)*x*[relation(n-1),0] and (n-1)*[0,0,relation(n-2)] and check the relative lengths of those; the two portions should be the same and so should be subtractable (and then divided by n)
T
T 2013-10-8
Yes they're the same...not sure what else to do.
Walter Roberson
Walter Roberson 2013-10-8
Could you confirm that even though the two portions are the same length, that when you ask to subtract them you are getting an error?
If
c = sum( (c1 - c2)/n )
then it follows that
c = (sum(c1) - sum(c2)) / n
and that form might not give the subtraction problem.
T
T 2013-10-8
Actually I avoided splitting:
c = ((2*n-1)*x*[relation(n-1),0] - (n-1)*[0,0,relation(n-2)])/n;
because using c1 and c2 does not generate the expected polynomial as you would expect here:
http://en.wikipedia.org/wiki/Legendre_polynomials
Walter Roberson
Walter Roberson 2013-10-8
What I was referring to as c1 was
(2*n-1)*x*[relation(n-1),0]
and for c2 it was
(n-1)*[0,0,relation(n-2)
breaking this into two assignments allows a breakpoint to be put in to confirm that the two are the same size before doing the subtraction
c = sum( (c1 - c2) / n )
which in turn is algebraically
c = (sum(c1) - sum(c2)) / n
you might need to simplify() the result for it to look sensible.
Is there a reason you are not using the explicit formulation involving the "choose" operator and x^k ?
T
T 2013-10-8
That's right. And for some reason using c1 and c2 will not generate the same polynomials for n where if you used
c = ((2*n-1)*x*[relation(n-1),0] - (n-1)*[0,0,relation(n-2)])/n;
The point was to do this using recursive programming.
T
T 2013-10-11
编辑：T 2013-10-11
Actually I don't think I need to do any of the above:
for n=0:inf
m=n+1;
int( ((2*n-1)*x*[relation(n-1),0] - (n-1)*[0,0,relation(n-2)])/n )* ((2*m-1)*x*[relation(m-1),0] - (m-1)*[0,0,relation(m-2)])/m),x,-1,1)
end
This will just give me zero. But is it possible to determine the round-off error? Or at least get the maximum value of n before the loop crashes?
I can set the recursion limit using:
set(0, 'RecursionLimit', 100000) but it keeps crashing still.

请先登录，再进行评论。

更多回答（1 个）

Matt J
Matt J 2013-10-6
编辑：Matt J 2013-10-6

0 个投票

Get rid of the "syms x". You are doing only numerical manipulations with the polynomial so there is no apparent need to have it in symbolic form, and that includes for plotting. Example
t=linspace(-1,1,1000);
p=[1 2 1]; %polynomial coefficients
plot(t,polyval(p,t))
Also, why use trapsum instead of the builtin trapz?

3 个评论
显示 1更早的评论 隐藏 1更早的评论

T
T 2013-10-6
I have removed syms x and ran with n = 8
However I still receive the error:
??? Undefined function or method 'relation' for input arguments of type 'double'.
How would I be able to plot the polynomial if the output for c returns an array with x as undefined?!
Matt J
Matt J 2013-10-6
It's the function relation() that can't be found. Check your path.
T
T 2013-10-6
Okay now I receive the error:
??? Undefined function or variable 'x'.
Error in ==> rlegendre at 9
c = ((2*n-1)*x*[rlegendre(n-1),0] - (n-1)*[0,0,rlegendre(n-2)])/n;
This happens when I remove syms x.

请先登录，再进行评论。

类别

在 帮助中心 和 File Exchange 中查找有关 Common Operations 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by