==============================================================================
SECTION A, QUESTION 1
=====================
f:=1/(1+x^2+sin(2*x)^2);
g:=diff(f,x);
fsolve(f+g=0,x=0.2);
fsolve(f+g=0,x=0.7);
fsolve(f+g=0,x=1.6);
fsolve(f+g=0,x=2.1);
fsolve( f:=int(subs(x=t,f),t=0..x), x=1 );
answers: 0.1085146471, 0.6472186476, 1.620476564, 2.065798378
0.6157566214
==============================================================================
SECTION A, QUESTION 2
=====================
function y=qua2a(x)
y=1./(1+x.^2+sin(2*x).^2);
function y=qua2b(x)
y=qua2a(x)-quad(@qua2a,0,x);
fzero(@qua2b,0.5);
gives answer 0.6158 in format short
0.615756650046297 in format long
if you want to improve the accuracy you need to change the function qua2b!!
==============================================================================
SECTION A, QUESTION 3
=====================
function z=qua3(x,y)
z=zeros(size(x));
for i=1:size(x,1)
for j=1:size(x,2)
if x(i,j)<1 & y<1+x(i,j) & y>-1-x(i,j)
z(i,j)=sqrt(x(i,j)^2+y^2);
end
end
end
dblquad(@qua3,-1,1,-2,2)
==============================================================================
SECTION A, QUESTION 4
=====================
function y=qua4a(x,p)
y=-(sin(x(1)+x(3))+sin(x(2)+x(3))+p*sin(x(1)-x(2)))/(1+x(1)^2+x(2)^2+x(3)^2);
function y=qua4b(p)
fminsearch(@qua4a,[0,0,0],[],p)
==============================================================================
SECTION A, QUESTION 5
=====================
with(LinearAlgebra):
a1:=<<6,0,0>|<0,1,0>|<0,0,-5>>;
a2:=<<1/2,-1/2,1>|<3,-8,1>|<1/2,-6,1>>;
M:=a2.a1.MatrixInverse(a2);
(if I now write Eigenvectors(M) I get the eigenvalues and eigenvectors in a
different order - in fact there is not a unique answer to this question!)
==============================================================================
SECTION B, QUESTION 1
=====================
a)
function z=qub1a(P)
r12=norm(P(1,:)-P(2,:));
r23=norm(P(2,:)-P(3,:));
r31=norm(P(3,:)-P(1,:));
z=sqrt((r12-r23)^2+(r23-r31)^2+(r31-r12)^2);
b)
function equts=qub1b(P,epsilon)
% epsilon is the error threshhold
n=size(P,1);
equts=[]; % this will hold the answer
% each row will give the 3 vertices of a triangle
for i=1:(n-2)
for j=1(i+1):(n-1)
for k=(j+1):n
z=qub1a([P(i,:);P(j,:);P(k,:)]);
if z<eps
equts=[equts;[P(i,:),P(j,:),P(k,:)]];
end if
end
end
end
==============================================================================
SECTION B, QUESTION 2
=====================
function y=qub2a(a,b,x)
% a is a vector of the form [1 0 a_{n-2} a_{n-1} ... a_0]
% b similar .....
% the quick way to do this:
ax=polyval(a,x)
bx=polyval(b,x)
y=abs(ax-bx)./sqrt(1+abs(ax).^2+abs(bx).^2);
% polyval(a,x) evaluates the polynomial at the points x (could be a vector)
% if you don't know polyval you have to write your own
% ax=a(1)*ones(size(x));
% for i=1:n
% ax=ax.*x+a(i+1);
% end
% bx=b(1)*ones(size(x));
% for i=1:n
% bx=bx.*x+b(i+1);
% end
Even though the question asked for the above function only to work for
scalar x, it works also for vector x.
Now for large |x|, the function in the integrand behaves like
|a(2)-b(2)|/\sqrt(2)x^2 , with an integral from L to infinity of
|a(2)-b(2)|/\sqrt(2)L. So you might try
L=1e5/|a(2)-b(2)|;
quad(@qub2a,-L,L)
better than this would be
L=1e5/abs(a(2)-b(2));
quad(@qub2a,-L,L) + sqrt(2)*1e-5
other ways to do it: integrate from 0 up to some L and then L to 2L - if the contribution
of the latter is small enough ignore it. Otherwise look at contribution from 2L to 3L etc.
L=10;
pos=quad(@qub2a,0,L);
add=pos;
i=1;
while add/pos>1e-5
add=quad(@qub2a,i*L,(i+1)*L);
pos=pos+add;
i=i+1
end
neg=quad(@qub2a,-L,0);
add=neg;
i=1;
while add/neg>1e-5
add=quad(@qub2a,-(i+1)*L,-i*L);
neg=neg+add;
i=i+1
end
==============================================================================
SECTION B, QUESTION 3
=====================
First of all - you MUST plot f,g to see where you should search for the minimum.
If they cross, there is nothing to do!!!
You must choose ranges for the plots until you get a good idea of what is
going on. Suppose we have done this and have found
(1) everything interesting happens in the range [a,b]
(2) x=x0 is a good starting approximation for locating both the minimum
vertical distance and minimal total distance.
with(plots);
zf:=plot(f,x=a..b);
zg:=plot(g,x=a..b);
For the minimal vertical distance:
xx:=fsolve( diff(f-g,x)=0, x=x0);
minvdist:=subs(x=xx,abs(f-g));
z1:=plot( [xx, subs(x=xx,f)+t*(subs(x=xx,g)-subs(x=xx,f)), t=0..1] )
For the minimal total distance:
dis:=sqrt((x1-x2)^2+(subs(x=x1,f)-subs(x=x2,g))^2);
xxx:=fsolve( {diff(dis,x1)=0,diff(dis,x2)=0}, {x1=x0,x2=x0} );
mindist:=subs(xxx,dis);
z2:=plot( [ subs(xxx,x1+t*(x2-x1)), subs(xxx,subs(x=x1,f)+t*(subs(x=x2,g)-subs(x=x1,f))),
t=0..1] );
display(zf,zg,z1,z2);
The horizontal distance is much worse! Need to do min(|x1-x2|) subject to the
constraint f(x1)=g(x2). Usual technique is Lagrange multipliers, but not easy
to implement in Maple.....
==============================================================================
SECTION B, QUESTION 4
=====================
The most important thing here is to get what each of the variables
represents.
Input: p1,p2 - 2 vectors specifying opposing points of the octahedron
Other variables:
v vector joining p1 to p2 (computed via v:=p2-p1)
c center of the octahedron (computed via c:=(p1+p2)/2)
r half the length of the vector v (computed via r:=Norm(v,2)/2)
e1,e2 two vectors of unit length perpendicular to v
(if v is of the form (0,0,a) then e1:=<1,0,0>; e2:=<0,1,0>;
else e1:=<v[2],-v[1],0>; e2:=CrossProduct(v,e1);
e1:=e1/Norm(e1,2); e2:=e2/Norm(e2,2);)
circ this is a circle, center c, radius r, in the plane spanned by e1,e2
the other 4 vertices of the octahedron must lie on this circle
t is the parameter of the circle
v1,v2,v3,v4 the other 4 vertices of the octahedron. they must lie on circ,
we take them at angles t=0,Pi/2,Pi,3Pi/2 but you could
take them at angles t=a,a+Pi/2,a+Pi,a+3Pi/2 for any a
d1,d2,d3,d4,d5,d6,d7,d8,d9,d10,d11,d12 plots of the 12 lines of the octahedron
verts a plot of the 6 vertices of the octahedron
The main point is that the other 4 vertices of the octahedron have to lie
on a circle (of suitable radius) in the plane perpendicular to the line
joining the two given vertices. Once you get this - understanding the
commands should be quite routine.
==============================================================================
Back to course homepage
Back to my teaching homepage
Back to my homepage