Saturday, August 27, 2016

Euler's Formula: A matLab tutorial

This tutorial was developed in MatLab for my daughter Penelope as a companion to chapter 4 of Robert Crease's The Great Equations, on Euler's Equation. Euler was a legendary calculator who spent his entire life doing calculations. Imagine if he had had MatLab!

Contents


Basic MatLab


Stuff entered into MatLab is shown in the grey boxes; MatLab's responses are shown below.

a value

1
ans =

     1

an expression

1+1
ans =

     2

assigning to a variable

x=1+1
x =

     2

recalling a variable

x
x =

     2

; after means don't show the results

x=1+1;

a complicated expression showing *, exponentiation, grouping

x=1+2*2+3^2+(5-3)/1.7
x =

   15.1765

calling a named function

factorial(4)
ans =

    24

a row vector

x=[4,1,7,1]
x =

     4     1     7     1

transpose operator ' turns a row into a column

x=[4,1,7,1]'
x =

     4
     1
     7
     1

: operator generates regularly spaced vectors

0:5
-3:.1:2
ans =

     0     1     2     3     4     5


ans =

  Columns 1 through 7

   -3.0000   -2.9000   -2.8000   -2.7000   -2.6000   -2.5000   -2.4000

  ...
  Columns 50 through 51

    1.9000    2.0000

Factorial can be applied to a vector

factorial(0:5) % returns a vector of factorials
ans =

     1     1     2     6    24   120

Euler's series for e


See also wikipedia

Note MatLab notation ./ to mean elementwise division of a vector. If you just use / it means something else.

n=10; % generate 10 terms in series
1./factorial(0:n) % returns terms in series
ans =

  Columns 1 through 7

    1.0000    1.0000    0.5000    0.1667    0.0417    0.0083    0.0014

  Columns 8 through 11

    0.0002    0.0000    0.0000    0.0000

sum up terms

sum(1./factorial(0:n))
ans =

    2.7183

sum is same as e

e=exp(1)
e =

    2.7183
 

which also happens to be 

Limit

nn=1:10000;
limit=(1+1./nn).^nn; % converges very slowly
limitStart=limit(1:20)
limitEnd=limit(end)
limitStart =

  Columns 1 through 7

    2.0000    2.2500    2.3704    2.4414    2.4883    2.5216    ...


limitEnd =

    2.7181

plot series and limit for e

clf
x=0:n;
plot(x,series,'*--'); hold on; plot(x,cumsum(series),'*--'); plot(xlim,exp(1)*[1,1],'y'); plot(nn(1:10),limit(1:10),'c*-'); hold off;
legend('terms','sum','e','limit (scaled)','location','SE'); title('Euler''s series for e'); shg

generate powers of x using a vector exponent

x=5; % pick a value of x
x.^(0:5)
ans =

           1           5          25         125         625        3125

Euler's series for e^x

series=x.^(0:n)./factorial(0:n)
series =

  Columns 1 through 7

    1.0000    5.0000   12.5000   20.8333   26.0417   26.0417   21.7014

  Columns 8 through 11

   15.5010    9.6881    5.3823    2.6911

sum it up, get exp(x) (almost -- need more terms)

sum(series)
exp(x) % e^x
clf; plot(0:n,series); hold on; plot(0:n,cumsum(series)); plot(xlim,exp(x)*[1,1],'y--'); hold off; shg
ylim(ylim*1.2);
title(sprintf('Euler series for e^{%3.2f}',x));
ans =

  146.3806


ans =

  148.4132

exponential and log

x=-5:.1:5;     subplot(1,2,1); plot(x,exp(x)); title('y=e^x=exp(x)');
x=[.1,1:100];  subplot(1,2,2); plot(x,log(x)); title('y=log_e(x)=ln(x)');

other bases

clf
x=-1:.1:2;     subplot(1,2,1); plot(x,[2.^x;exp(x);10.^x]); legend('2^x','e^x','10^x'); title('Exponentials with different bases');
x=[.1,1:100]; subplot(1,2,2); plot(x,[log2(x);log(x);log10(x)]'); legend('log_2(x)','log_e(x)','log_{10}(x)','location','SE');  title('Logs with different bases');
shg

sin and cos

degreesPerRadian=360/(2*pi); % multiplier to convert radians to degrees
radiansPerDegree=1/degreesPerRadian;
subplot(2,2,1);
x=-pi:.1:2*pi; % in radians
plot(x*degreesPerRadian,sin(x)); title('y=sin(x)'); xlabel('degrees'); settick('X',90*[-2:4]);
line(xlim,[0,0],'color','k');  line([0,0],[-1,1],'color','k');
% line(90*[1,1],[-1,1],'color','y'); line(180*[1,1],[-1,1],'color','y'); line(270*[1,1],[-1,1],'color','y');
settick('X',90*[-2:4]);
subplot(2,2,2);
plot(x*degreesPerRadian,cos(x)); title('y=cos(x)'); xlabel('degrees'); settick('X',90*[-2:4]);
line(xlim,[0,0],'color','k');  line([0,0],[-1,1],'color','k');

subplot(2,2,3); plot(x*degreesPerRadian,[sin(x);cos(x)]'); legend('sin(x)','cos(x)'); xlabel('degrees');  settick('X',90*[-2:4]); shg
line(xlim,[0,0],'color','k');  line([0,0],[-1,1],'color','k');
title('sin(x) = cos(x - 90^o)');

subplot(2,2,4); plot(x.^2.*sin(x.^2)); title('Random complicated expression: x^2sin(x^2)');
line(xlim,[0,0],'color','k');  line([0,0],[-1,1],'color','k');

sin vs cos

close all; w=-10:.1:10;
plot(sin(w),cos(w)); axis square; xlabel('sin(w)'); ylabel('cos(w)');  shg

tan

clf
x=-360:360;
plot(x,tan(x/degreesPerRadian)); ylim([-10,10]); shg

animate circle, sin, cos, tan

if(false) % don't use with MatLab's publish; won't work
    close all; w=-10:.1:10;% fullscreen;
    clear M;
    sinw=sin(w); cosw=cos(w);
    x=0:360; xr=x*radiansPerDegree; sinx=sin(xr); cosx=cos(xr); tanx=tan(xr); shg
    for i=0:360;
        r=i*radiansPerDegree; sinr=sin(r); cosr=cos(r); tanr=tan(r);
        subplot(2,2,1);  plot(sinw,cosw,'-',[0,cosr],[0,sinr],'-',cos(r),sinr,'*'); axis square; xlabel('sin(w)'); ylabel('cos(w)');
        subplot(2,2,2);  plot(x,sinx,'-',i,sin(r),'*'); title('sin(w)');
        subplot(2,2,3);  plot(x,cosx,'-',i,cos(r),'*'); title('cos(w)');
        subplot(2,2,4);  plot(x,tanx,'-',i,tanr,'*'); ylim([-5,5]); title('tan(w)');
        drawnow;
        M(i+1)=getframe(gcf);
    end
    % create an movie of the animation
    vw= VideoWriter('pers/Trig.avi');
    open(vw);
    writeVideo(vw,M);
    close(vw) % then use MakeAGif.com to convert to animated gif
end
irrationals & complex numbers

sqrt(-1)
ans =

   0.0000 + 1.0000i

MatLab defines symbol i to mean sqrt(-1)

clear i; % restore i to its mathematical meaning rather than a variable
i
ans =

   0.0000 + 1.0000i

i^2 = -1

i^2
ans =

    -1

you can make any complex # you want

3*i+2
ans =

   2.0000 + 3.0000i

One way to make alternating signs

(-1).^(0:n)
ans =

     1    -1     1    -1     1    -1     1    -1     1    -1     1

Generate an odd numbered series

1:2:10
ans =

     1     3     5     7     9

Euler series for sin

n=10;
ii=1:2:n; % odd term indices
x=.8; % a random # to find the sin of
sign=(-1).^(0:length(ii)-1)
num=x.^ii
denom=factorial(ii)
series=sign.*num./denom
sumSeries=sum(series)
sinx=sin(x)
clf; plot(ii,series,'*--',ii,cumsum(series),'*--',[min(ii),max(ii)],sinx*[1,1],'y');
legend('terms','sum','sin(x0','location','SE'); title(sprintf('Euler''s series for sin(%3.2f)',x)); shg
sign =

     1    -1     1    -1     1


num =

    0.8000    0.5120    0.3277    0.2097    0.1342


denom =

           1           6         120        5040      362880


series =

    0.8000   -0.0853    0.0027   -0.0000    0.0000


sumSeries =

    0.7174


sinx =

    0.7174

Euler series for cos

n=10;
ii=0:2:n;  % even term indices
x=.8; % a random # to find the sin of
sign=(-1).^(0:length(ii)-1)
num=x.^ii
denom=factorial(ii)
series=sign.*num./denom
sumSeries=sum(series)
cosx=cos(x)
sign =

     1    -1     1    -1     1    -1


num =

    1.0000    0.6400    0.4096    0.2621    0.1678    0.1074


denom =

           1           2          24         720       40320     3628800


series =

    1.0000   -0.3200    0.0171   -0.0004    0.0000   -0.0000


sumSeries =

    0.6967


cosx =

    0.6967

Euler series for e^ix

x=.736; % a random x
n=10
ii=0:n
num=(i*x).^ii
denom=factorial(ii)
series=num./denom
sumSeries=sum(series)
cosx=cos(x)
sinx=sin(x)
clf; plot(ii,real(series),'*--',ii,imag(series),'*--',ii,real(cumsum(series)),'^--',ii,imag(cumsum(series)),'^--');
legend('real component of terms','imaginary component  of terms','real component of sum','imaginary component  of sum','location','SE');
title(sprintf('Euler''s series for e^{%3.2fi}',x)); shg
n =

    10


ii =

     0     1     2     3     4     5     6     7     8     9    10


num =

  Columns 1 through 4

   1.0000 + 0.0000i   0.0000 + 0.7360i  -0.5417 + 0.0000i   0.0000 - 0.3987i

  Columns 5 through 8

   0.2934 + 0.0000i   0.0000 + 0.2160i  -0.1590 + 0.0000i   0.0000 - 0.1170i

  Columns 9 through 11

   0.0861 + 0.0000i   0.0000 + 0.0634i  -0.0466 + 0.0000i


denom =

  Columns 1 through 6

           1           1           2           6          24         120

  Columns 7 through 11

         720        5040       40320      362880     3628800


series =

  Columns 1 through 4

   1.0000 + 0.0000i   0.0000 + 0.7360i  -0.2708 + 0.0000i   0.0000 - 0.0664i

  Columns 5 through 8

   0.0122 + 0.0000i   0.0000 + 0.0018i  -0.0002 + 0.0000i   0.0000 - 0.0000i

  Columns 9 through 11

   0.0000 + 0.0000i   0.0000 + 0.0000i  -0.0000 + 0.0000i


sumSeries =

   0.7412 + 0.6713i


cosx =

    0.7412


sinx =

    0.6713

Euler series for e^i*pi

x=pi % a not-so-random x
n=20
ii=0:n
num=(i*x).^ii
denom=factorial(ii)
series=num./denom
sumSeries=sum(series)
cosx=cos(x)
sinx=sin(x)
fprintf('QED!!!!!');
x =

    3.1416


n =

    20


ii =

  Columns 1 through 13

     0     1     2     3     4     5     6     7     8     9    10    11    12

  Columns 14 through 21

    13    14    15    16    17    18    19    20


num =

   1.0e+09 *

  Columns 1 through 4

   0.0000 + 0.0000i   0.0000 + 0.0000i  -0.0000 + 0.0000i   0.0000 - 0.0000i

  Columns 5 through 8

   0.0000 + 0.0000i   0.0000 + 0.0000i  -0.0000 + 0.0000i   0.0000 - 0.0000i

  Columns 9 through 12

   0.0000 + 0.0000i   0.0000 + 0.0000i  -0.0001 + 0.0000i   0.0000 - 0.0003i

  Columns 13 through 16

   0.0009 + 0.0000i   0.0000 + 0.0029i  -0.0091 + 0.0000i   0.0000 - 0.0287i

  Columns 17 through 20

   0.0900 + 0.0000i   0.0000 + 0.2828i  -0.8886 + 0.0000i   0.0000 - 2.7916i

  Column 21

   8.7700 + 0.0000i


denom =

   1.0e+18 *

  Columns 1 through 7

    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000

  Columns 8 through 14

    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000

  Columns 15 through 21

    0.0000    0.0000    0.0000    0.0004    0.0064    0.1216    2.4329


series =

  Columns 1 through 4

   1.0000 + 0.0000i   0.0000 + 3.1416i  -4.9348 + 0.0000i   0.0000 - 5.1677i

  Columns 5 through 8

   4.0587 + 0.0000i   0.0000 + 2.5502i  -1.3353 + 0.0000i   0.0000 - 0.5993i

  Columns 9 through 12

   0.2353 + 0.0000i   0.0000 + 0.0821i  -0.0258 + 0.0000i   0.0000 - 0.0074i

  Columns 13 through 16

   0.0019 + 0.0000i   0.0000 + 0.0005i  -0.0001 + 0.0000i   0.0000 - 0.0000i

  Columns 17 through 20

   0.0000 + 0.0000i   0.0000 + 0.0000i  -0.0000 + 0.0000i   0.0000 - 0.0000i

  Column 21

   0.0000 + 0.0000i


sumSeries =

  -1.0000 - 0.0000i


cosx =

    -1


sinx =

   1.2246e-16

i.e., exp(i*pi)+1=0 QED!!!


Series for ln(x)=log_e(x)

See bottom formula of first table here
x=1.357;
y=(x-1)/(x+1)
ii=1:2:n % odd terms only
series=2*y.^ii./ii
sumSeries=sum(series)
log(x)
y =

    0.1515


ii =

     1     3     5     7     9


series =

    0.3029    0.0023    0.0000    0.0000    0.0000


sumSeries =

    0.3053


ans =

    0.3053

Leibnitz' series for pi

See wikipedia
n=10000 % this one converges verrrrry slowly
ii=1:2:n; display(ii(1:10)); % odd terms only. show only first few.
sign=(-1).^(0:length(ii)-1); display(sign(1:10));
series=4*sign./ii; display(series(1:10));
sumSeries=sum(series)
pi
k=20;
plot(ii(1:k),series(1:k),'-',ii(1:k),cumsum(series(1:k)),'-',[min(ii),k],pi*[1,1],'y');
legend('terms','sum','pi'); shg
n =

       10000


ans =

     1     3     5     7     9    11    13    15    17    19


ans =

     1    -1     1    -1     1    -1     1    -1     1    -1


ans =

  Columns 1 through 7

    4.0000   -1.3333    0.8000   -0.5714    0.4444   -0.3636    0.3077

  Columns 8 through 10

   -0.2667    0.2353   -0.2105


sumSeries =

    3.1414


ans =

    3.1416

No comments:

Post a Comment