======================================================================
Part A, Question 1
==================
f:=sin(x)-x+x^3/1000;
g:=diff(f,x);
Digits:=5;
The commands below work...but it is not easy to read off the ranges
directly from the graph, it requires a little playing around....
fsolve( g=0, x=-1..1 ) gives 0
fsolve( g=0, x=5..6 ) 5.8278
fsolve( g=0, x=6..7 ) 6.8176
fsolve( g=0, x=11..12 ) 11.632
fsolve( g=0, x=13..14 ) 13.683
fsolve( g=0, x=15..20 ) 17.373
fsolve( g=0, x=20..22 ) 20.711
fsolve( g=0, x=22..25 ) 22.944
======================================================================
Part A, Question 2
==================
Make a file q2.m with the following content:
function z=q2(x,y,p)
z=zeros(size(x));
for i=1:size(x,1)
for j=1:size(x,2)
if abs(y)<abs(x(i,j)) & x(i,j)^2+y^2<1
z(i,j)=sin(sqrt(1+x(i,j)^2+p*y^2));
end
end
end
and then type dblquad(@q2,0,1,-1,1,[],[],p)
(where p is the value you want)
======================================================================
Part A, Question 3
==================
with(student) loads the "student" library.
I1:=... defined I1 to be the integral shown, which is
NOT evaluated (as Int is used and not int).
I2:=... executes the substitution sin(x)=u in I1.
the integral is now with respect to the variable u
I3:=.. integrates I2 by parts, where u^2 is differentiated.
"simplify" simplifies the final result
I4:=.. integrates I3 by parts, where u is differentiated.
"simplify" simplifies the final result
value(I1) asks Maple to find the result of integrating
I1 directly (i.e. what we would have got from "int"
instead of "Int")
clearly it can be seen to be the same as I4.
The computation being done here is a comparison of the
result of the "step-by-step" evaluation of the original
integral I1 with the immediate answer.
======================================================================
Part A, Question 4
==================
Make a file q4.m with the content
function z=q4(x)
z=2*x(1)^2+x(2)^2-x(1)/(1+x(1)^2+2*x(2)^2)-1;
and then do [a b]=fminsearch(@q4,[0,0])
(get answer a=[0.2172 0.0000] b=-1.1131)
to make graph make a new file q4.m with the content
function z=q4(x,y)
z=2*x.^2+y.^2-x./(1+x.^2+2*x.^2)-1;
and then do the commands
x=[-2:0.05:2];
X=ones(81,1)*x;
Y=X';
mesh( X,Y,q4(X,Y) )
see graph.
======================================================================
Part A, Question 5
==================
f:=proc(a,N)
return( a/N*sum(1/(a+i^2), i=1..N ));
end proc;
plot3d( f(a,N), a=1..4, N=10..100 );
f as we have defined it here is perfectly well defined for integer N
and noninteger N (in which case it will compute the sum up to floor(N)).
Whenever you plot a function you always sample it at a discrete set
of points, so you should not worry that the original f was defined
only for integer N.
see graph.
======================================================================
Part B, Question 1
==================
We will assume that S takes the form of a SET of LISTS
( [1,2,3], [3,4,5] } or similar
a)
q11:=proc(S)
local d,dnew,dims,i,j,k;
d:=0;
dims:=nops(op(1,S));
for i from 1 to nops(S) do
for j from (i+1) to nops(S) do
dnew:=0;
for k from 1 to dims do
dnew:=dnew+(op(k,op(i,S))-op(k,op(j,S)))^2;
end do;
dnew:=evalf(sqrt(dnew));
if dnew>d then
d:=dnew;
end if;
end do;
end do;
return(d);
end proc;
b)
q12:=proc(S)
local i,x,y,z;
x:=0;
y:=0;
z:=0;
for i from 1 to nops(S) do
x:=x+op(1,op(i,S)):
y:=y+op(2,op(i,S)):
if nops(op(i,S))=3 then
z:=z+op(3,op(i,S)):
end if;
end do:
if nops(op(1,S))=2 then
return( [x/nops(S),y/nops(S)] );
else
return( [x/nops(S),y/nops(S),z/nops(S)] );
end if;
end proc;
c)
q13:=proc(S)
local d,dnew,dims,i,M;
M:=q12(S);
d:=0;
dims:=nops(op(1,S));
for i from 1 to nops(S) do
dnew:=0;
for k from 1 to dims do
dnew:=dnew+(op(k,op(i,S))-op(k,M))^2;
end do;
dnew:=evalf(sqrt(dnew));
if dnew>d then
d:=dnew;
end if;
end do;
return(d);
end proc;
d)
q14:=proc(S)
local D,R;
D:=q11(S);
R:=q13(S);
if D>2*R then
return("diameter bigger than twice the radius");
else
return("diameter not bigger than twice the radius");
end if;
end proc;
======================================================================
Part B, Question 2
==================
a)
make a file q21.m:
function y=q21(x,p)
y=x./(3+p(1)*sin(x)+p(2)*cos(2*x));
make a file q22.m:
function z=q22(p)
z=quad(@q21,0,2*pi,[],[],p);
and now type [a b] = fminsearch( @q22, [0,0] )
get a = [-1.0170 -0.0211 ]
b = 6.2236
Suppose in the integral we make errors of order e.
We cannot expect the value of the function at the minimum to be
more accurate than this, and thus the error in the values of p and
q (in the vector a) will be of the order sqrt(e). If you use a
tolerance e in the computation of the integral, do not make the
"xtol" for fimnsearch more than sqrt(e) !!
b)
make a file q2b1.m:
function y=q2b1(x,p,alpha)
y=x.^alpha./(3+p(1)*sin(x)+p(2)*cos(2*x));
make a file q2b2.m:
function z=q2b2(p,alpha)
z=quad(@q2b1,0,2*pi,[],[],p,alpha);
make a file q2b3.m
function z=q2b3(alpha)
z=zeros(size(alpha));
for i=1:size(alpha,1)
for j=1:size(alpha,2)
[q z(i,j)]=fminsearch(@q2b2,[0,0],[],alpha(i,j));
end
end
and now type quad(@q2b3,1,2)
get 12.8672 (though the accuracy is probably limited)
======================================================================
Part B, Question 3
==================
a)
At the maximum value of x for which there is a solution to the
equation the graph is vertical, i.e. dy/dx = infinity.
Remember that if f(x,y) is constant then dy/dx=-(df/dx)/(df/dy).
So we need df/dy=0. This gives 2y-x=0, so y=x/2 at this point.
f:=y^2-x*y+x^3;
f1:=subs(y=x/2,f);
fsolve(f1=2,x); gives an answer about 1.349
To find where there are changes in the number of y's we need to
find the critical points, i.e. where dy/dx=0 => df/dx=0.
Since df/dx=-y+3x^2 we
f2:=subs(y=3*x^2,f);
xs:=[fsolve(f2=2,x)];
this gives x vals of crit points, -0.637, 0.750
to get the relevant y values do
y1:=fsolve(subs(x=op(1,xs),f)=2)
gives -1.85, 1.22 , the latter is the relevant solution
y2:=fsolve(subs(x=op(2,xs),f)=2)
gives -0.936, 1.69, the latter is the relevant solution
b)
Method 1:
There are changes in the number of x's at points where dy/dx=infinity,
i.e. df/dy=0. So we use
f:= ... (write in the function) ;
fsolve( {f=0,diff(f,y)=0}, {x,y} );
There are changes in the number of y's at points where dy/dx=0
i.e. df/dx=0. So we use
f:= ... (write in the function) ;
fsolve( {f=0,diff(f,x)=0}, {x,y} );
BUT FSOLVE ONLY GIVES A SINGLE SOLUTION FOR A SYSTEM....so you will
have to do this many times over, using the graph as a guide how to start
fsolve each time.
Method 2:
You could use a for loop to find how many solutions there are for
each x and for each y:
f:= ... (write in the function) ;
for x from -10 to 10 by 0.1 do
x;
fsolve( f, y);
end do;
similarly switching y and x (but do restart first!)
======================================================================
Part B, Question 4
==================
a)
function z=q41(x,p,q)
% find p(x)
valpx=0;
for i=1:size(p,1)
valpx=valpx+p(i)*x^(size(p,1)-i);
end
% find q(x)
valqx=0;
for i=1:size(q,1)
valqx=valqx+q(i)*x^(size(q,1)-i);
end
% find roots
z=roots([1 valpx valqx]');
% eliminate complex roots
if imag(z(1))~=0
z=[];
end
b)
function q42(p,q,a,b)
h=(b-a)/100;
n=1; % number of points on graph so far
for i=0:100
x=a+i*h;
ys=q41(x,p,q);
for j=1:size(ys,1)
X(n)=x;
Y(n)=ys(j);
n=n+1;
end
end
plot(X,Y,'b.')
c) Assume the coefficients of the polynomials are given in a
matrix with d columns which we will call P. The function q42.m
hardly needs to be changed, just replace "p,q" in two places by
"P". The modification of the function q41.m will be much easier
if we take the coefficients in P to be in the reverse order
from usual. Then we can write:
function z=q41(x,P)
for i=1:size(P,2)
p(i)=0;
for j=1:size(P,1)
p(i)=p(i)+P(j,i)*x^(j-1);
end
end
% find roots
y=roots([1; p']);
%remove complex roots
z=[];
for i=1:size(y,1)
if imag(y(i))==0
z=[z;y(i)];
end
end
[ it really works! try
q42([-1 0 0 0; -2 0 0 3; -3 2 4 1]',-3,3)
and compare against
implicitplot( y^3-y^2+(-2+x^3)*y+(-3+2*x+4*x^2+x^3)=0,x=-3..3,y=-10..10)
in Maple. ]
======================================================================
Back to course homepage
Back to my main teaching page
Back to my main page