Using a special version of the Maple Kernel, the Toolbox performs symbolic analysis using character strings, rather than numerical analysis on numerical arrays.
maple('sum(''1/k!'', ''k''=0..infinity)') %note the use of double quotesThe Symbolic Toolbox is a collection of tools for manipulating and solving symbolic expressions. When you pass a command to Matlab to perform a symbolic operation, it asks Maple to do it and return the result to the Matlab command window.
Using the Matlab Toolbox we can write
H=sym(15,15,'1/(i+j-1)'); x=sym(15,1,'1'); b=symmul(H,x); flops(0); linsolve(H,b); flopsOn the other hand, in Matlab syntax we have
H=hilb(15);x=ones(15,1); b=H*x; flops(0); x=H\b; flops cond(H)
What do you learn from a direct comparison of exact and floating point solutions of this linear system?
H=hilb(3) mlr=eig(H)
Now we redo the calculation using Maple
H=maple('hilbert','3') r=eigensys(H)From a practical point of view this answer is useless, since we don't really have any idea what is the size of these numbers. Lets convert to floating point.
mar=numeric(r) maple('evalf',r,50)
Notice this last answer, while correct to any realistic digits, is still rather strange since H is symmetric and we would expect to be returned real numbers like Matlab did. The point is that sometimes floating point answers are better than exact or floating point approximations of exact answers.
Here are some examples:
diff('cos(x)') % differentiate cos(x) with respect to x M=sym('[a,b;c,d]') % create a symbolic matrix M determ(M) % find the determinant of the matrix M
'cos (x)'or
'1/(1+x^2)'
diff('cos(t)')in simple cases the program can guess the variable to differentiate with respect to without being told.
int('1/(1+x^2)')
f='x^n'or
f='sin(a*t+b)'then more general syntax for diff and int are
diff(f,'x')or
diff(f,'n')Note the difference
diff(f,'x',k)where k is an integer differentiates f, k times.
int(f,'t')finds the indefinite integral.
int(f,'t',a,b)finds the definite integral from a to b.
For example try
int('2*x^2*(19+12*x^2)/(7*(x^2+1))','x',0,1) numericThe first ans gives a better answer from a mathematical point of view.
Here are some exercises:
x^2*sin(x^3)*e^{2x}
x^6*log(cos(x^2)*x)
x^3*e^{x^4}
(ax+b)/(cx+d)
atan(x^2+1)
x^7
1/x
log(x)*sqrt(x)
x^7 % on (0,1)
1/x % on (1,2)
sqrt(x)*log(x) % on (0,1)
exp (-x^2) % on (0,infinity)
f='exp(i*theta)' diff(f,'theta',2)gives
-exp(i*theta)
f='1/(5+4*cos(x))' ezplot(f) f2=diff(f,2) ezplot(f2) axis([-2*pi 2*pi -2 2])here axis is a matlab command to set an axis range. Here is another example
f3=diff(f2) pretty(f3)Lets try to simplify this
f4=simple(f3) pretty(f4)the command simple tries to simplify a complex expression. It will try many different Maple options. Lets look at root finding in the toolbox:
z=solve(f3)the command solve finds zeros of the function. The answer z is a symbolic matrix. To find numeric values issue the next command.
x=numeric(z)numeric without an argument numerically evaluates the last expression. Sometimes you might want to pick off, say, the second answer in z. The next line shows one way to do this
s=sym(z,2,1)sym is the symbol manipulator. Note that z is a 3 by 1 matrix so we type the pair 2,1 to indicate that entry. Often you want to substitute a value into a symbolic expression. This can be done in several ways. One way is as follows:
subs(f3,s)Now lets try to recover f
g=int(int(f2))Note this does not look like f. To check our work lets take f-g and see what we get. Now f and g are strings in matlab with different lengths, so you cannot simply take f-g. To carry out symbolic difference (symsub for symbolic subtraction), you would type the following:
d=symsub(f,g)Note the answer is not zero but, rather, as we know from calculus the answer can differ by a constant.
simple(ans)You should use mhelp to look at the following related topics symadd, symmul, symdiv, sympow, symop, symrat, symvar, symsize.
f='(x+y)^10*(2*u+3*v)^7*(x-2*y)^3' g=expand(f)the command expand multiplies all the terms out. As simple examples, try expanding
(x-1)*(x-2)*(x-3) x*(x*(x-6)+11)-6 a*(x+y) exp(a+b) n*(n-1)!
Now consider the reverse process factoring:
factor(g)
Try factoring
x^3-6*x^2+11x-6 x^6+1 x^4+4 x^4-2*x^2+1
for n=1:9 p=symadd(sympow('x',n),1); disp([p '= ' factor(p)]) end
for n=1:17 N=setstr('1'*ones(1,n)); fn=factor(N); disp([N,' = ',fn]) end
simplify(f)here f must be a maple expression. Try letting f be the following expressions and use simplify
(1-x^2)/(1-x) (1/a^3+6/z^2+12/a+8)^{1/3} e^xe^y cos^2(x)+sin^2(x)
2*cos^2(x)-sin^2(x) cos^2(x)-sin^2(x) cos(x)+i*sin(x) cos(3\arccos(x))
n=10 f1='2*sin(j*pi*x)*x^2*(1-x)^2'; for j=1:n f=subs(f1,j,'j'); ss=['x=0','..','1']; vvv=(maple('int',f,ss)); y(j)=eval(vvv) end
clear n=6; v=ones(1,n); t=sym(v); f='cos(k*acos(x))'; for k=1:n t=sym(t,1,k,subs(f,k,'k')) endNote that these don't look like polynomials. Lets fix that using simple and expand
tcheb=sym(ones(1,n)); for j=1:n %tcheb=sym(tcheb,1,j,simple(sym(t,1,j))) tcheb=sym(tcheb,1,j,expand(sym(t,1,j))) endThe Tchebshev polynomials can also be built by recursion. Lets do this using matlab polynomials which, you may recall are given by vectors, i.e., if . The recusion formula for the Tchebyshev polynomilas is
Note that the polynomial is given in matlab by [1, 2, 3]. In order to express the polynomial xp(x) we would write [1, 2, 3,0]. Using this information we have
t0=1; t1=[1 0]; for j=1:(n-1) eval(['t',int2str(j+1),' = [2*t',int2str(j),' 0]-[0 0 t',int2str(j-1),']']) endIf you really like to see the x's we can a symbolic toolbox feature that turns Matlab polynomials into regular variable polynomials.
for j=1:n eval(['tchebp',int2str(j),' = poly2sym(t',int2str(j),')']) end
Now lets plot our answers: First lets plot the matlab polynomials t1, t2, t3 t4 t5 t6
clear y xvec=0:.05:1; for j=1:6 eval(['y(:,',int2str(j),')=polyval(t',int2str(j),',xvec)'';']) end plot(xvec,y) axis([0 1 -1.5 1.5]);
Before ending this topic we should mention that many important special functions including the standard orthogonal polynomials are built into Maple but they are only available with the Extended symbolic Toolbox.
The Characteristic equation for this problem can be written as
Consider the case and . Note that, using x for this equation can also be written in terms of sin's and cos's. Plotting the resulting function
ezplot('(1*2-x^2)*(sin(x)/x)+(1+2)*cos(x)',[.0001,4*pi])we easily see that there are infinitely many zeros interlacing with multiples of .
Lets find the first 10 zeros.
k0=1 k1=2 f=subs('(k0*k1-x^2)*(sin(x)/x)+(k0+k1)*cos(x)',k0,'k0'); f1=subs(f,k1,'k1'); for j=1:10 ss=['x=Pi*',int2str(j-1),'..',int2str(j),'*Pi']; vvv=(maple('fsolve',f1,ss)); mu(j)=eval(vvv); end mu'
Consider the one dimensional heat equation on the interval 0 to 1:
The solution to this problem can be written in terms of a Dirichlet series:
where
Imagine that we can sample the temperature at a single point on the rod at a sequence of times to obtain output . Using this information we want to approximate the Fourier coefficients and hence approximate the initial condition . Defining the matrices
we arrive at a system of equations:
Assuming that is irrational, then the matrix P is invertible and, as a generalized Vandermonde matrix, A is invertible. Therefore, in principle, we can solve for X and obtain:
With this we obtain an approximation to f, namely,
Lets consider an example:
In a physical experiment we would measure the outputs using some type of sensing device but in the present case we will simply approximate the solution using the Dirichlet series formula for the solution. To do this we must compute the Fourier coefficients of the initial data. Lets do it using the symbolic toolbox and using matlab:
First lets use Maple
nn=20 for j=1:nn f=subs('2*x^2*(1-x)^2*sin(j*pi*x)',j,'j'); ss=['x=0','..','1']; vvv=(maple('int',f,ss)); cma(j)=eval(vvv); end
We could also have typed
maple('g:=(j)->int(2*x^2*(1-x)^2*sin(j*pi*x),x=0..1)') sg=maple('[seq(g(j),j=1..20)]') cma=numeric(sg) % these are the Fourier coefficients of f
Now lets use Matlab:
We can build a function file SPMquot
coefn.m", for example,
function y=coefn(x) global jj y=2*x.^2.*(1-x).^2.*sin(jj*pi*x);
Now write the program SPMquot
fcml"
global jj nn=20 % use Tchebyshev quadrature nc=50; % the number of Tchebyshev nodes xc=cos((2*(1:nc)-1)*pi/(2*nc)); %nodes for Tchebyshev quadrature % you can use any interval with (a , b) for the integration, here we use a=0; b=1; for j=1:nn jj=j; cml(j)= ((b-a)*pi)/(2*nc)* sqrt(1-xc.^2)*coefn((b+a+(b-a)*xc)/2)'; end cml' % these are the Fourier coefficients of fNotice that the result was almost instantaneous.
Lets compare the results
format long [cma' cml']
Check the Fouier coefficients
nn=8; x=0:.01:1; J=(1:nn); cmly=cml(1:nn)*sin(pi*J'*x); plot(x,cmly,x,(x.*(1-x)).^2)Looks like a good approximation
Lets look at numerical solution of the heat problem on
n=20; t=[0:.05:2]'; x=0:.05:1; J=(1:n); v=cml'; % build a column vector of Fourier coefs vv=v(:,ones(size(x))); %this is a neat trick to build v*ones(size(x)) Y=exp(-pi^2*t*J.^2)*(vv.*sin(J'*pi*x)); surf(x,t,Y)That's pretty easy. Okay back to solution of our inverse problem
Now suppose we sample the temperature at n=8 time values, say at the times , at the x value .
x_0=1/pi; nn=20; n=8; T=(1:n).^2/(n^2); J=(1:nn); Y=exp(-pi^2*T'*J.^2)*(cml(1:nn).*sin(J*pi*x_0))'
n=8; x_0=1/pi; V=1:n; A=exp(-pi^2*(T)'*(V.^2)); P=diag(sin(V*pi*x_0)); AP=A*P; X=AP\Y
nn=8; x=0:.01:1; J=(1:nn); y=X'*sin(pi*J'*x); plot(x,y,x,(x.*(1-x)).^2)Looks bad, lets try a few less coefficients
nn=3; x=0:.01:1; J=(1:nn); y=X(1:nn)'*sin(pi*J'*x); plot(x,y,x,(x.*(1-x)).^2)Well that looks a bit more like .
Now lets give Maple a chance: We could rebuild the above matrices using the Maple toolbox, for example,
maple('f:=(i,j)->exp(-pi^2*i^2*j^2/8^2)') SA=maple('matrix(8,20,f)') vsg=maple('matrix(20,1,[seq(g(j)*sin(j),j=1..20)])') symmul(SA,vsg)
We find that results cannot be to long.
A simpler way is as follows
SA=sym(A); SP=sym(P); SY=sym(Y); XS=symmul(inverse(SP),linsolve(SA,SY)) %SX=maple('evalf',XS) SX=numeric(XS)
Lets try our plot with nn=3 as above
nn=3; x=0:.01:1; J=(1:nn); sy=SX(1:nn)'*sin(pi*J'*x); plot(x,sy,x,(x.*(1-x)).^2)Lets try nn=2
nn=2; x=0:.01:1; J=(1:nn); sy=SX(1:nn)'*sin(pi*J'*x); plot(x,sy,x,(x.*(1-x)).^2)Well it might be okay for government work but not really very good.
Once again the exact answer is of no use and an attempt to obtain a floating point approximation to it fails. The theory of ill-posed inverse problems would suggest that we ``regularize'' the problem and solve as follows:
n=8; x_0=1/pi; V=1:n; A=exp(-pi^2*(T)'*(V.^2)); P=diag(sin(V*pi*x_0)); RX=inv(P)*((A'*A +.000065*eye(n))\(A'*Y)); [cml(1:n)' X SX RX]
nn=8; x=0:.01:1; J=(1:nn); sy=RX(1:nn)'*sin(pi*J'*x); plot(x,sy,x,(x.*(1-x)).^2)
The pictures better but at least we can use all the approximate Fourier coefficients in our approximation.
To be fair in this problem please redo the whole thing using the following sample times
n=8; T=(1:n).^2/(pi^2*n^2);and
T=(1:n).^2/(n^4);
You will see that for the last T you can use all eight approximate coefficients. But the regularized solution does very bad. This is because you must adjust the parameter used to multiply times the identity (in this case change .000065 to zero.
f='sinh(a*sqrt(sqrt(-1)*x))/sinh(sqrt(sqrt(-1)*x))'on . We are really only interested in the real part which is oscillatory but decays like . From results in Fourier Analysis it can be shown that the real part of the integral of f from zero to infinity is zero.
First, lets use Matlab's quad8 to approximate the integral. To this end build two function files ``yr.m"
function y=yr(x) global aa z=x+eps; % we put in eps to take care of the 0/0 when x=0 y=real(sinh(aa*sqrt(sqrt(-1)*z))./sinh(sqrt(sqrt(-1)*z)));and ``yi.m''
function y=yi(x) global aa z=x+eps; y=imag(sinh(aa*sqrt(sqrt(-1)*z))./sinh(sqrt(sqrt(-1)*z)));Now build an m-file
SPMquot
int_y.m"
global aa aa=input(' global variable aa in (0,1) = '); b(1)=0; int_yr(1)=0; int_yi(1)=0; disp(' b int_yr int_yi') disp('--------------------------------------------') for j=2:7 b(j)=10^j; int_yr(j)=quad8('yr',b(j-1),b(j))+int_yr(j-1); int_yi(j)=quad8('yi',b(j-1),b(j))+int_yi(j-1); fprintf('%8.0f \t %8.8f\t %8.8f\n',b(j),int_yr(j),int_yi(j)) endReturn to Matlab and type
SPMquot
int_y" return. Try ``aa=.5'' and ``aa=.9''
Well we might guess that the integral of the real part is zero and the imag part is minus two hundred when ``aa=.9''.
Now lets try the symbolic toolbox: First the real part
maple('evalf(Re(Int(2*sinh(.9*i^.5*v)/sinh(i^.5*v)*v,v=0..10000)),8)')and now the imaginary part
maple('evalf(Im(Int(2*sinh(.9*i^.5*v)/sinh(i^.5*v)*v,v=0..10000)),8)')
Now the fact is that the integral
can be expressed in terms of the Zeta function by
So it is clear that for every ``aa'' the real part is zero and we can compute the result using the symbolic toolbox to obtain
maple('-i*(Zeta(0,2,1/2*(1-.9))-Zeta(0,2,1/2*(1+.9)))/2')