Om Sri Guru Raghavendra
HW1
HW2
HW3
HW4
HW1-Updated
HW2-Updated
HW5
HW5-Updated
HW6
HW7
<br\>
Mtg1
(1) Plotting a given function ![{\displaystyle f(x)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/202945cce41ecebb6f643f31d119c514bec7a074)
[edit | edit source]
Ref: Lecture Notes p.2-2
(i) Find
(ii)Plot
(i) By L'Hospital's rule,
|
|
|
|
(ii) MATLAB Code to generate the plot
clc;
clear all;
x=0:0.1:1;
for i=1:length(x)
if x(i)==0
y(i)=1;
else
y(i)=(exp(x(i)) - 1)/ x(i);
end
end
plot(x,y,'r*-');
title('Plot of (exp(x)-1)/x vs x');
xlabel('x');
ylabel('f(x)');
Plot:
vs
(7) Expressing a function
using Taylor Series
[edit | edit source]
Ref: Lecture Notes p.6-1
(i)Construct a Taylor series of
around
for ![{\displaystyle \displaystyle n=0,1,2......,10.}](https://wikimedia.org/api/rest_v1/media/math/render/svg/54912abf8c0fab37e63c015e2e97f1dae4e69f63)
(ii)Plot the series for each
(iii)Estimate the maximum
at
(i)
order polynomial is given by,
(a)
for
and
|
|
|
(b)
for
and
|
|
|
|
(c)
for
and
|
|
|
|
(d)
for
and
|
|
|
|
(e)
for
and
|
|
|
|
(f)
for
and
|
|
|
|
(g)
for
and
|
|
|
|
(h)
for
and
|
|
|
|
(i)
for
and
|
|
|
|
(j)
for
and
|
|
|
|
(k)
for
and
|
|
|
|
(ii) MATLAB Code for generating the plots for :
clc;
clear all;
syms x;
f = sin (x);
x_0 = pi/4;
y=0:0.01:pi;
n=0:1:10;
for j=1:length(n)
for i=1:length(y)
if n(j)== 0
z(i,j)=((y(i)-x_0)^n(j) / (factorial(n(j))))* subs(diff(f,x,n(j)),x,x_0);
else
z(i,j)=z(i,j-1) + ((y(i)-x_0)^n(j) / (factorial(n(j))))* subs(diff(f,x,n(j)),x,x_0);
end
end
end
for k=1:5
if n(k)==0
plot (y,z(:,1),'r-');
hold on;
end
hold on;
if n(k)==1
plot (y,z(:,2),'b-');
hold on;
end
if n(k)==2
plot (y,z(:,3),'g-');
legend('n=2');
end
if n(k)==3
plot (y,z(:,4),'k-');
hold on;
end
if n(k)==4
plot (y,z(:,5),'m-');
hold on;
end
if n(k)==5
plot (y,z(:,6),'c-');
hold on;
end
legend('n=0','n=1','n=2','n=3','n=4','n=5');
end
title('sin (x) vs x for n=0,1,2,3,4 and 5');
xlabel('x');
ylabel('sin (x)');
hold on;
for k=6:length(n)
if n(k)==6
plot (y,z(:,7),'r-');
hold on;
end
hold on;
if n(k)==7
plot (y,z(:,8),'b-');
hold on;
end
if n(k)==8
plot (y,z(:,9),'g-');
legend('n=2');
end
if n(k)==9
plot (y,z(:,10),'k-');
hold on;
end
if n(k)==10
plot (y,z(:,11),'m-');
hold on;
end
legend('n=6','n=7','n=8','n=9','n=10');
end
title('sin (x) vs x for n=6,7,8,9,and 10');
xlabel('x');
ylabel('sin (x)');
hold on;
Plot:
for
for
(iii)Error Estimation
![{\displaystyle \displaystyle R_{n+1}(x)={\frac {(x-x_{0})^{n+1}}{(n+1)!}}f^{n+1}(\xi )}](https://wikimedia.org/api/rest_v1/media/math/render/svg/0c0f372d3ce2fa69d2fb33d5f18a0fbb5f089d7c)
Since
for all
(a)
for
,
and
|
|
|
|
(b)
for
,
and
|
|
|
|
(c)
for
,
and
|
|
|
|
(d)
for
,
and
|
|
|
|
(e)
for
,
and
|
|
|
|
(f)
for
,
and
|
|
|
|
(g)
for
,
and
|
|
|
|
(h)
for
,
and
|
|
|
|
(i)
for
,
and
|
|
|
|
(j)
for
,
and
|
|
|
|
(k)
for
,
and
|
|
|
|
--Subramanian Annamalai 23:58, 26 January 2010 (UTC)
1. Run Kessler's code to reproduce the table of constants and the values of polynomials at t=1 obtained by Kessler.
2. Explain Kessler's code by giving comments for every line.
3. Obtain
by understanding Kessler's code line by line.
Ref: Lecture Notes p.30-2 and p.37-1
Ref: Trapezoidal rule error
Part I : Generation of the Table
To generate the table given in Kessler's paper, the following modifications
are to be carried out.
1. The first half of the code is the main code.
Give that a name called kessler.m
2. The second half of the code is a subroutine with the name fracsum.
This should not be part of the file kessler.m
We need to create a new file called fracsum.m in the same folder
that contains the file kessler.m
3. Set n=8. This is to ensure that we get the constants
and polynomials
.
4. The most important part is that the subroutine fracsum should end with the line,
dsum=dsum/div;
5. The following lines have been appended to kessler.m to generate the table.
format long;
q_c = (cn./cd); % Will generate the quotients of c
q_p = (pn./pd); % Will generate the quotients of p
%Printing results
disp(c); disp(q_c);
disp(p); disp(q_p);
Result
The constants
Constant |
Numerator |
Denominator |
Quotient
|
![{\displaystyle c_{1}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/77b7dc6d279091d354e0b90889b463bfa7eb7247) |
-1 |
1 |
-1.00000000000000
|
![{\displaystyle c_{3}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/b1dc52bfbaf6e577fbed72a716068f4533700bd3) |
1 |
6 |
0.16666666666667
|
![{\displaystyle c_{5}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/407ede56ab53d2b4d945a82cfde0146e896189fe) |
-7 |
360 |
-0.01944444444444
|
![{\displaystyle c_{7}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/67c5bf4c3a3a4d9883e64a807be8f898def72c29) |
31 |
15120 |
0.00205026455026
|
![{\displaystyle c_{9}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/ed85b9f062a13f13f8a2e1e4835a074cf7ea2236) |
-127 |
604800 |
-0.00020998677249
|
![{\displaystyle c_{11}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/c28e5cf59854c4b6d0c8ffdc6a964b495b0ee770) |
73 |
3421440 |
0.00002133604564
|
![{\displaystyle c_{13}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/8259df723480fab04dde715947b5030f1536235c) |
-1414477 |
653837184000 |
-0.00000216334744
|
![{\displaystyle c_{15}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/144e4ab3def0874d2c872de0ec4f80364abbda2d) |
8191 |
37362124800 |
0.00000021923271
|
![{\displaystyle c_{17}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/115f44221cc848718155d957f01bbc409c8ee596) |
-16931177 |
762187345920000 |
-0.00000002221393
|
The polynomials evaluated at t=1
Polynomial |
Numerator |
Denominator |
Quotient
|
![{\displaystyle p_{2}(1)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/ba9f9f21439e0ae32bd40251962d2d7c3463739d) |
-1 |
3 |
-0.33333333333333
|
![{\displaystyle p_{4}(1)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/1ab185ae537e86adff385d242b26e8489569e6b6) |
1 |
45 |
0.02222222222222
|
![{\displaystyle p_{6}(1)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/2d35a2e4261b8f8775b9cccdf823bfbf5ee66957) |
-2 |
945 |
-0.00211640211640
|
![{\displaystyle p_{8}(1)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/05938e665a9f1181dcf96e5ab06c9e68937686d7) |
1 |
4725 |
0.00021164021164
|
![{\displaystyle p_{10}(1)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/2d26281bd8d9cff358193ca929d1bb43b1c5aa05) |
-2 |
93555 |
-0.00002137779916
|
![{\displaystyle p_{12}(1)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/a1b997267b5e999652e9573f8509ea34fdc95913) |
1382 |
638512875 |
0.00000216440428
|
![{\displaystyle p_{14}(1)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/05158f8475434a617e4dda348c55eaa14350f646) |
-4 |
18243225 |
-0.00000021925948
|
![{\displaystyle p_{16}(1)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/27171f6d1fabc1eee6a1a298dd40d52895412da9) |
3617 |
162820783125 |
0.00000002221461
|
Part II: Line by Line Explanation of Kessler's MATLAB Code - Adding Comments
function [c,p]=traperror(n)
%traperror is a user-defined function which takes in 'n' as the input argument
%and gives the constants and the coefficients of t in the polynomial p.
f=1; g=2;
cn=-1; cd=1;
%cn is the numerator of the constants in the polynomial
%cd is the denominator of the constants in the polynomial
%cn is initially set to 1 and cd to -1 because c1 = -1, i.e,(-1/1).
for k=1:n
%The process is repeated for as many number of polynomials as desired.
f=f.*g.*(g+1);
% If we observe carefully we find that the denominators of the constants c1, c2, c3... are odd-numbered.
% f generates these odd-numbered factorials.
% For example, 1! =1 , 3! =6 and so on.
[newcn,newcd]=fracsum(-1*cn,cd.*f);
%The function fracsum is called to get the values of the numerator
%and denominator of the constants.
%Again on observing that, once if we find the constant c1,
% the constant c3 can be found from the equation, c1/1! + c3/3! =0
% Then using c1 and c3, the constant c5 can be found using,
% c1/5! + c3/3! + c1/1! = 0 and so on.
% Hence the negative of the constant and the corresponding factorials are needed.
% These are in fact the arguments of the function fracsum.
cn=[cn;newcn]; cd=[cd;newcd];
% The current returned by the function fracsum and the
% previous values of cn and cd are stored
f=[f;1];
g=[2+g(1);g];
%The next value of f and g are generated.
% The denominator factorials of the polynomials evaluated at t=1 differ from
% the previous values of the constants only by 1.
% Hence f and g are so adjusted to evaluate p(1) for all n.
[newpn,newpd]=fracsum((g-1).*cn,f.*cd);
%The function fracsum() returns the new values of the numerator and denominator
%of the coefficient of t which is now stored in newpn and newpd respectively.
pn(k,1)=newpn; pd(k,1)=newpd;
% Here we assign p(1) of the numerator and denominator for varying n.
end
c=int64([cn cd]);
%This gives the numerator and denominator of the constants.
p=int64([pn pd]);
%This gives the numerator and denominator of the coefficients of t.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [nsum,dsum]=fracsum(n,d)
%fracsum is a user-defined function that takes in all the
%previous values of the constants and generates the new values.
div=gcd(round(n),round(d));
% The input arguments n and d are rounded of to the nearest integer values
% This is carried out using the in-built "round" function
%Then the greatest common divisor of n and d is found using the in-built 'gcd' function.
% This is done in order that we have the constants in fractions - Only then can we call it as cn and cd.
% Otherwise we will end up in a decimal.
n=round(n./div);
% The numerator is divided by the GCD to end up in a whole number
d=round(d./div);
% The numerator is divided by the GCD to end up in a whole number
dsum=1;
%The denominator of the constant is initialized to 1.
for k=1:length(d)
dsum=lcm(dsum,d(k));
%The least common factor of the new values of the denominator
%of the constant and the previous value is calculated.
%Here since we have to add all the fractions LCM is calculated.
end
%Then dsum gives the new value of the denominator of the constant.
nsum=dsum*sum(n./d);
div=gcd(round(nsum),round(dsum));
% Again the gratest common divisor is found between the numerator and denominator
% of the constants using the 'gcd' and 'round' functions as before.
% As before this is to ensure that we end up in a whole number and not a fraction.
nsum=nsum/div;
%nsum gives the new value of the numerator of the constant.
dsum=dsum/div;
%dsum gives the new value of the denominator of the constant.
|
Part III : Generating
<br\>
|
![{\displaystyle \displaystyle (Eq.1)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/53bc710cf3779357c030fecd81c7bf815aae4396)
|
|
![{\displaystyle \displaystyle (Eq.2)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/549b9b71fb3dd781fd437c65c44e3b0bff3821db)
|
|
![{\displaystyle \displaystyle (Eq.3)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/67029a61096aa9304e425bf080fbeedf36390330)
|
|
![{\displaystyle \displaystyle (Eq.4)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/698bb882efb26bc61f16ec2c4a36578d59e90561)
|
|
![{\displaystyle \displaystyle (Eq.5)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/285fe6908ded39d6692cb8e1660f47e457751d81)
|
|
![{\displaystyle \displaystyle (Eq.6)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/e93f27477d6127264f6858b0c7bc4758905a298d)
|
|
![{\displaystyle \displaystyle (Eq.7)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/0ebf085a76a7258f58a5a9d18e9797b4fbf80b75)
|
<br\>
Hence we need to find out the constants,
|
<br\>
So we start off with setting
And then we need to have n=3 ,because,
|
|
|
|
<br\>
Consider the first iteration: i.e, k=1
|
yields
f= 1*2*3 =6
|
Hence,
|
|
|
<br\>
|
yields
fracsum(-1*1, 1*6)
i.e, fracsum(-1,6)
|
Hence,we have
|
|
|
<br\>
Now looking into the fracsum(n,d) function,
![{\displaystyle \displaystyle div=gcd(round(n),round(d))}](https://wikimedia.org/api/rest_v1/media/math/render/svg/45ce1f2dda4a54347fd6efbdf7f5de4a93926535)
|
will yield div = 1.
because,
round(-1) = -1
round(6) = 6
and hence,
gcd(-1,6) = 1.
|
<br\>
Now,
![{\displaystyle n=round(n./div);}](https://wikimedia.org/api/rest_v1/media/math/render/svg/3e39a7ee53ec970a360a37c545f615173f53fd83)
will give n=1
becasue round(1/1) =1
|
and,
![{\displaystyle \displaystyle d=round(d./div);}](https://wikimedia.org/api/rest_v1/media/math/render/svg/112b40a70ff63ac417ad93b392c70d31918d96fc)
will give d=6
becasue round(6/1) =6
|
Thus,
|
|
|
Since length(d) = 1
and
lcm (1,6) = 6,
we have, dsum =6.
|
Thus,
|
|
|
![{\displaystyle \displaystyle div=gcd(round(nsum),round(dsum));}](https://wikimedia.org/api/rest_v1/media/math/render/svg/afb0d348d7a2ead2df784478ebfdc44e0d8d60e9)
|
will yield,
nsum = 6* sum(1/6) = 1
and
|
<br\>
and,
div = 1
since,
round(nsum) = round(1) = 1
round(dsum) = round(6) = 6
gcd(1,6) =1
|
Now,
|
will yield ,
Now going back,
|
Thus,
|
![{\displaystyle \displaystyle \therefore c_{3}={\frac {1}{6}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/cc0d3ecb52ae2c987f735a195267626ee8537a20)
|
|
<br\>
Next, we have,
|
Hence,
cn = [-1;1]; cd=[1;6]
|
This is followed by,
;
which means, g = [2+2, 2] = [4;2]
|
Now consider the second iteration - i.e, k=2
Consider the second iteration: i.e, k=2
|
yields
f= (6*4*5; 1*2*3)
|
Hence,
|
|
|
<br\>
|
yields
fracsum(-1*-1, 1*120)
i.e, fracsum(1,120)
|
and,
yields
fracsum(-1*1, 6*6)
i.e, fracsum(-1,36)
|
Hence,we have
|
|
|
<br\>
Now looking into the fracsum(1,120) function,
![{\displaystyle \displaystyle div=gcd(round(n),round(d))}](https://wikimedia.org/api/rest_v1/media/math/render/svg/45ce1f2dda4a54347fd6efbdf7f5de4a93926535)
|
will yield div = 1.
because,
round(1) = 1
round(120) = 120
and hence,
gcd(1,120) = 1.
|
<br\>
Now,
![{\displaystyle n=round(n./div);}](https://wikimedia.org/api/rest_v1/media/math/render/svg/3e39a7ee53ec970a360a37c545f615173f53fd83)
will give n=1
becasue round(1/1) =1
|
and,
![{\displaystyle \displaystyle d=round(d./div);}](https://wikimedia.org/api/rest_v1/media/math/render/svg/112b40a70ff63ac417ad93b392c70d31918d96fc)
will give d=120
becasue round(36/1) =120
|
and,
<br\>
Now looking into the fracsum(-1,36) function,
![{\displaystyle \displaystyle div=gcd(round(n),round(d))}](https://wikimedia.org/api/rest_v1/media/math/render/svg/45ce1f2dda4a54347fd6efbdf7f5de4a93926535)
|
will yield div = 1.
because,
round(-1) = -1
round(36) = 36
and hence,
gcd(-1,36) = 1.
|
<br\>
Now,
![{\displaystyle n=round(n./div);}](https://wikimedia.org/api/rest_v1/media/math/render/svg/3e39a7ee53ec970a360a37c545f615173f53fd83)
will give n=-1
becasue round(-1/1) =-1
|
and,
![{\displaystyle \displaystyle d=round(d./div);}](https://wikimedia.org/api/rest_v1/media/math/render/svg/112b40a70ff63ac417ad93b392c70d31918d96fc)
will give d=6
becasue round(36/1) =36
|
Thus,
|
|
|
Since length(d) = 2
and
For d(1),
dsum = lcm (1,120) = 120,
and once again for d(2),
dsum = lcm (120,36) = 360
Finally we have,
dsum = 360.
|
Thus,
|
|
|
![{\displaystyle \displaystyle div=gcd(round(nsum),round(dsum));}](https://wikimedia.org/api/rest_v1/media/math/render/svg/afb0d348d7a2ead2df784478ebfdc44e0d8d60e9)
|
will yield,
nsum = 360* sum(1/120 - 1/36) = -7
and
|
<br\>
and,
div = 1
since,
round(nsum) = round(-7) = -7
round(dsum) = round(360) = 360
gcd(-7,360) =1
|
Now,
|
will yield ,
Now going back,
|
Thus,
|
![{\displaystyle \displaystyle \therefore c_{5}=-{\frac {7}{360}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/93828be5d0655b10cda033e4520a9d0b34a527a8)
|
|
<br\>
Next, we have,
|
Hence,
cn = [-1;1;7]; cd=[1;6;360]
|
This is followed by,
;
which means, g = [2+4;4; 2] = [6;4;2]
|
<br\>
Again following the same procedure for Iteration 3, i.e for k=3
we will obtain,
|
![{\displaystyle \displaystyle \therefore c_{7}={\frac {31}{15120}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/c860fd58575b2b7bd645fd4e3c3fa43e31c2b23f)
|
|
<br\>
<br\>
On substituting, the values of
|
in (Eq.1) through (Eq.7) we obtain,
|
|
|
|
|
|
|
|
<br\>