Underdamped harmonic oscillator

(1)

DE, StudentNumber: 91

Underdamped harmonic oscillator

Let

a= 0.4, b= 1.04, t₀= 0, t₁= 6, (1)

ic0 = 2, ic1=−1, (2)

y⁰⁰(t) +ay⁰(t) +by(t) =f(t). (3)

The first order form of (3) is d dt

y v

=A y

v

+B(f(t)) =

0 1

−b −a y v

+

0 1

(f(t)), (4)

wherey⁰=v. The initial condition (2) is

~ y(0) =

y(0) v(0)

= ic₀

ic1

. (5)

Mathematica. {a,b,t0,t1,ic0,ic1}={?,?, ?,?, ?,?}

de={y’’[t]+a y’[t]+b y[t]==f[t]}

ic={y[t0]==ic0,y’[t0]==ic1}

A={{0,1},{-b,-a}}

Octave. clear all;

pkg load symbolic;

global a=? b=? A=[0,1;-b, -a];

t0=?; t1=?; ic0=?; ic1=?;

syms y(t);

de=diff(y,t,2)+a*diff(y,t)+b*y(t);

function ydot = velocity (yvec,t) global a b A;

ydot=A*yvec;

endfunction N=100;

1, A

Exercise. Let f(t) = 0, (2) be true. How much isy(t₁) ?

Mathematica. desol=NDSolve[Join[de/.{f[t]->0},ic],y,{t,t0,t1}]

ans=(y[t1]/.desol)[[1]]

ans=Replace[y[t1],desol][[1]]

Octave. ts = linspace (t0, t1, N);

ys = lsode ("velocity", [ic0; ic1], ts);

disp ("The value of y(t1) is:"), disp (ys(N,1)) Answers. A: 0.628891

Correct answer. A

(2)

2, Aa

Exercise. Letf(t) = 0, (2) be true. LetU_t₁_,t₀(~z) =~y(t₁), where~y(t) is the solution of (4) with the initial condition

~

y(t0) =~z. How much is the first component of

U_t₁_,t₀( ic₀

ic₁

).

disp ("The value of (yvec(t1))(1) is:"), disp (ys(N,1)) Answers. A: 0.628891

Correct answer. A 3, Ab

Exercise. Letf(t) = 0, (2) be true. LetUt1,t0(~z) =~y(t1), where~y(t) is the solution of (4) with the initial condition

~

y(t₀) =~z. Compute U_t₁_,t₀(~z), if~z=n₁(0.2,0)^T +n₂(0,0.2)^T, n_1,2∈ {0,1,2,3,4,5}. Plot these points!

Mathematica. m=N[MatrixExp[t1*A]]

Graphics[Join[{PointSize[Large],Point[{0,0}]},{PointSize[Medium]}, Table[Point[{x,y}],{x,0,1,0.2},{y,0,1,0.2}]//Flatten,{Red},

Table[Point[m.{x,y}],{x,0,1,0.2},{y,0,1,0.2}]//Flatten]]

Octave. figure

np=6; dx=1/(np-1); black=zeros(2,np^2); red=zeros(2,np^2);

u=expm(t1*A);

for iy=1:np for jx=1:np

black(:,1+(iy-1)*np+(jx-1))=[(jx-1)*dx;(iy-1)*dx];

red(:,1+(iy-1)*np+(jx-1))=u*[(jx-1)*dx;(iy-1)*dx];

endfor endfor;

hold on;

scatter(red(1,:),red(2,:),"r");

scatter(black(1,:),black(2,:),"b");

hold off;

title("[y(t0),y’(t0)] -> [y(t1),y’(t1)]")

(3)

Answers.

Correct answer. A 4, Ac

Exercise. f(t) = 0, (2) is true. Solve the following algebro-differential equation!

y⁰=v, v⁰ =−spring−f riction, spring=by, f riction=av.

How much isspring(t1) ?

Mathematica. ade={y’[t]==v[t],v’[t]==-friction[t]-spring[t], friction[t]==a v[t],spring[t]==b y[t],y[0]==ic0,v[0]==ic1};

adesol=NDSolve[ade,{y,v,spring,friction},{t,t0,t1}];

(spring[t1]/.adesol)[[1]]

Octave. x0=[ ic0; ic1; b*ic0; a*ic1 ];

xdot0=[ ic1; -b*ic0-a*ic1; b*ic0; a*ic1 ];

function res = resudials (x,xdot,t) global a b;

y=x(1); v=x(2); s=x(3); f=x(4);

ydot=xdot(1); vdot=xdot(2); sdot=xdot(3); fdot=xdot(4);

res(1)=ydot-v;

res(2)=vdot+s+f;

res(3)=s-b*y;

res(4)=f-a*v;

endfunction N=100;

ts = linspace(t0,t1,N);

[x, xdot] = daspk("resudials", x0, xdot0, ts);

disp ("The value of spring(t1) is:"); disp (x(N,3));

(4)

Answers. A: 0.654046 Correct answer. A 5, B

Exercise. f(t) = 0. The general solution of (3) is

y(t) =C1e^λ¹^t+C2e^λ²^t, =(λ₁)>0.

How much is=(C₁), if (2) is satisfied?

Mathematica. de2=Join[de/.{f[t]->0},ic]

sol2=DSolve[de2,y,t]

expSol=TrigToExp[sol2[[1,1,2,2]]]

Octave. rewrite( dsolve(de==0,y(t0)==ic0,diff(y,t)(t0)==ic1), ’exp’) Answers. A: 0.3

Correct answer. A 6, Ba

~

y(t) =C1e^λ¹^t~v1+C2e^λ²^t~v2, =(λ₁)>0, ~v1 = (1, s21)^T, ~v2 = (1, s22)^T How much is=(C₁), if (10) is satisfied?

Mathematica. esys=Eigensystem[A]//N evals=esys[[1]]

evects=esys[[2]]

If[Im[evals[[1]]]>0,v1=evects[[1]];v2=evects[[2]], v1=evects[[2]];v2=evects[[1]]]

icond={ic0,ic1}

vn1=v1/v1[[1]]

vn2=v2/v2[[1]]

csol=Solve[Table[(C1 vn1+C2 vn2)[[i]]==icond[[i]],{i,2}],{C1,C2}]

Im[C1/.csol]

Octave. [evects, evals] = eig (A);

vn1=[1;evects(2,1)/evects(1,1)];

S=[vn1 vn2];

inv(S)*[2;-1];

(inv(S)*[2;-1])(1);

disp ("Im(C1) is:");

disp( imag((inv(S)*[2;-1])(1) ));

Answers. A: 0.3 Correct answer. A 7, C

y(t) =e^αt(C₁cosωt+C₂sinωt).

How much isC₁, if (2) is satisfied?

(5)

(y[0]/.sol3)[[1]]

Octave. dsolve(de==0,y(t0)==ic0,diff(y,t)(t0)==ic1)

Answers. A: 2.

Correct answer. A 8, D

Exercise. f(t) = 0, (2), (3) true. Plot y(t) ! Mathematica. de4=Join[de/.{f[t]->0},ic]

sol4=NDSolve[de4,y,{t,t0,t1}]

Plot[Evaluate[y[t]/.sol4],{t,t0,t1}]

figure

plot(ts,ys(:,1));

title ("y(t)");

Answers.

1 2 3 4 5 6

-1.0 -0.5 0.5 1.0 1.5 2.0

Correct answer. A 9, E

Exercise. f(t) = 0, (2), (3) is true. Plot y(t), y⁰(t)-t on a single figure!

Plot[Evaluate[{y[t],y’[t]}/.sol5],{t,t0,t1}]

figure

plot(ts,ys(:,1));

title ("y(t),y’(t)");

(6)

Answers.

1 2 3 4 5 6

-1 1 2

Correct answer. A 10, F

Exercise. f(t) = 0, (2), (3) are true. Plot the

γ : [t₀, t₁]→R², γ(t) = (y(t), v(t))^T parametric curve!

ParametricPlot[Evaluate[{y[t],y’[t]}/.sol6],{t,t0,t1}]

figure

plot(ys(:,1),ys(:,2));

title ("t->[y(t),y’(t)]");

(7)

Answers.

-1.0 -0.5 0.5 1.0 1.5 2.0

-1.5 -1.0 -0.5 0.5 1.0

Correct answer. A 11, G

γ : [t0, t1]→R³, γ(t) = (t, y(t), v(t))^T parametric curve!

ParametricPlot3D[Evaluate[{t,y[t],y’[t]}/.sol7],{t,t0,t1}]

figure

plot3(ts,ys(:,1),ys(:,2));

title ("t->[t,y(t),y’(t)]");

(8)

Answers.

0

2

4

6 -1

0 1

2 -1

0

Correct answer. A 12, H

Exercise. Plot the (y, v)^T →A(y, v)^T vector field!

Mathematica. VectorPlot[A.{y,v},{y,-.2,.2},{v,-1,1}, VectorPoints->15,VectorStyle->Arrowheads[0.02]]

Octave. [hor,vert]=meshgrid(linspace(-0.3,0.3,10),linspace(-1,1,10));

figure

quiver (hor, vert, A(1,1)*hor+A(1,2)*vert,A(2,1)*hor+A(2,2)*vert);

title ("[y;v] -> A[y;v]");

(9)

Answers. ^-^1.0 ^-^0.5 ^0.0 ^0.5 ^1.0

-1.0 -0.5 0.0 0.5 1.0

Correct answer. A 13, Ha

Exercise. f(t) = 0. Plot the (y, v)^T →A(y, v)^T vector field and a solution of (4)!

Mathematica. sol=NDSolve[Join[de/.{f[t]->0},ic],y,{t,t0,3*t1}]

cPlot=ParametricPlot[Evaluate[{y[t],y’[t]}/.sol],{t,t0,3*t1}]

vPlot=VectorPlot[A.{y,v},{y,-1,1},{v,-1,1}, VectorPoints->8,VectorStyle->Arrowheads[0.04]]

Show[vPlot,cPlot]

Octave. figure hold on;

[hor,vert]=meshgrid(linspace(-1,1,10),linspace(-1,1,10));

plot(ys(:,1),ys(:,2));

title ("[y;v] -> A[y;v], t->[t,y(t),y’(t)]");

hold off;

(10)

Answers. ^-^1.0 ^-^0.5 ^0.0 ^0.5 ^1.0

-1.0 -0.5 0.0 0.5 1.0

Correct answer. A 14, I

Exercise. Compute U = exp(1.4A) ! How much isU11 ? Mathematica. ans9=MatrixExp[1.4 A]

ans9[[1,1]]

Octave. disp ("exp(1.4 * A) is:"), disp ( expm(1.4*A) ) Answers. A: 0.277416

Correct answer. A 15, J

Exercise. Compute U(t) = exp(tA)-t and plot its first column’s functions!

Mathematica. U=MatrixExp[t A]

Plot[{U[[1,1]],U[[2,1]]},{t,t0,t1}]

Octave. Uts=zeros(N,2,2);

for i=1:N

Uts(i,:,:)=expm(ts(i)*A);

endfor;

figure

plot(ts,Uts(:,1,1),ts,Uts(:,2,1)) title("first column of exp(tA)")

(11)

Answers.

1 2 3 4 5 6

-0.5 0.5 1.0

Correct answer. A 16, K

Exercise. f(t) = 0, (2) true, y(0) = 1, y⁰(0) = 0. Plot (y(t), y⁰(t))^T ! Mathematica. de11=Join[de/.{f[t]->0},ic]

Octave. ys = lsode ("velocity", [1; 0], ts);

figure

plot(ts,ys);

title ("y(t),y’(t), y(0)=1,y’(0)=0");

Answers.

1 2 3 4 5 6

-0.5 0.5 1.0

Correct answer. A 17, L

Exercise. Compute U(t) = exp(tA)-t and plot its second column’s functions!

Plot[{U[[1,2]],U[[2,2]]},{t,t0,t1}]

(12)

for i=1:N

endfor;

figure

plot(ts,Uts(:,1,2),ts,Uts(:,2,2)) title("second column of exp(tA)")

Answers.

1 2 3 4 5 6

-0.5 0.5 1.0

Correct answer. A 18, M

figure

plot(ts,ys);

title ("y(t),y’(t), y(0)=0,y’(0)=1");

Answers.

1 2 3 4 5 6

-0.5 0.5 1.0

Correct answer. A

(13)

Overdamped harmonic oscillator

Let

a= 0, b= 2, t0= 0, t1 = 2.5, (6)

ic₀= 0.8, ic₁ = 1.1, (7)

y⁰⁰(t) +ay⁰(t) +by(t) =f(t). (8)

y v

=A y

v

+B(f(t)) =

0 1

−b −a y v

+

0 1

(f(t)), (9)

~ y(0) =

y(0) v(0)

= ic₀

ic1

. (10)

Mathematica. {a,b,t0,t1,ic0,ic1}={?,?, ?,?, ?,?}

de={y’’[t]+a y’[t]+b y[t]==f[t]}

ic={y[t0]==ic0,y’[t0]==ic1}

A={{0,1},{-b,-a}}

Octave. clear all;

%graphics_toolkit ("gnuplot");

pkg load symbolic;

global a=? b=? A=[0,1;-b, -a];

t0=?; t1=?; ic0=?; ic1=-?;

syms y(t);

ydot=A*yvec;

endfunction N=100;

19, A

Exercise. Let f(t) = 0, (2) be true. How much isy(t1) ?

Correct answer. A 20, Aa

~

U_t₁_,t₀( ic0

ic₁

).

(14)

Exercise. Let f(t) = 0, (2) be true. Let U_t₁_/6,t₀(~z) = ~y(t1/6), where ~y(t) is the solution of (4) with the initial condition~y(t0) =~z. ComputeU_t₁_/6,t₀(~z), if~z=n1(0.2,0)^T+n2(0,0.2)^T, n1,2∈ {0,1,2,3,4,5}. Plot these points!

Mathematica. m=N[MatrixExp[(t1/6)*A]]

Octave. figure

u=expm((t1/6)*A);

endfor endfor;

hold on;

hold off;

title("[y(t0),y’(t0)] -> [y(t1),y’(t1)]")

(15)

Answers.

How much isspring(t₁) ?

y=x(1); v=x(2); s=x(3); f=x(4);

res(1)=ydot-v;

res(2)=vdot+s+f;

(16)

res(3)=s-b*y;

res(4)=f-a*v;

endfunction N=100;

y(t) =C1e^λ¹^t+C2e^λ²^t, λ1 > λ2. How much isC1, if (2) is satisfied?

sol2[[1,1,2,2]]

Octave. dsolve(de==0,y(t0)==ic0,diff(y,t)(t0)==ic1) Answers. A: 2.7

Correct answer. A 24, Ba

~

y(t) =C1e^λ¹^t~v1+C2e^λ²^t~v2, λ1> λ2, ~v1 = (1, s21)^T, ~v2 = (1, s22)^T How much is=(C₁), if (10) is satisfied?

Mathematica. esys=Eigensystem[A]//N evals=esys[[1]]

evects=esys[[2]]

If[evals[[1]]>evals[[2]], v1=evects[[1]];v2=evects[[2]], v1=evects[[2]];v2=evects[[1]]]

icond={ic0,ic1}

vn1=v1/v1[[1]]

vn2=v2/v2[[1]]

csol=Solve[Table[(C1 vn1+C2 vn2)[[i]]==icond[[i]],{i,2}],{C1,C2}]

C1/.csol

Octave. [evects, evals] = eig (A);

S=[vn1 vn2];

inv(S)*[2;-1];

(inv(S)*[2;-1])(1);

disp ("Im(C1) is:");

disp( imag((inv(S)*[2;-1])(1) ));

Answers. A: 2.7 Correct answer. A

(17)

25, D

figure

plot(ts,ys(:,1));

title ("y(t)");

Answers. ^0.5 ^1.0 ^1.5 ^2.0 ^2.5

0.4 0.6 0.8 1.0

figure

plot(ts,ys(:,1));

(18)

Answers.

0.5 1.0 1.5 2.0 2.5

-0.5 0.5 1.0

Correct answer. A 27, F

figure

plot(ys(:,1),ys(:,2));

title ("t->[y(t),y’(t)]");

(19)

Answers.

0.4 0.6 0.8 1.0

-0.4 -0.2 0.2

γ : [t₀, t₁]→R³, γ(t) = (t, y(t), v(t))^T parametric curve!

figure

plot3(ts,ys(:,1),ys(:,2));

title ("t->[t,y(t),y’(t)]");

(20)

Answers.

0

1

2

0.2 0.4

0.6 0.8

-0.4 -0.2

0.0 0.2

figure

title ("[y;v] -> A[y;v]");

(21)

Answers. ^-^1.0 ^-^0.5 ^0.0 ^0.5 ^1.0

-1.0 -0.5 0.0 0.5 1.0

Show[vPlot,cPlot]

plot(ys(:,1),ys(:,2));

hold off;

(22)

Answers. ^-^1.0 ^-^0.5 ^0.0 ^0.5 ^1.0

-1.0 -0.5 0.0 0.5 1.0

ans9[[1,1]]

Plot[{U[[1,1]],U[[2,1]]},{t,t0,t1}]

for i=1:N

endfor;

figure

(23)

Answers.

0.5 1.0 1.5 2.0 2.5

-0.5 0.5 1.0

figure

plot(ts,ys);

title ("y(t),y’(t), y(0)=1,y’(0)=0");

Answers.

0.5 1.0 1.5 2.0 2.5

-0.5 0.5 1.0

Plot[{U[[1,2]],U[[2,2]]},{t,t0,t1}]

(24)

for i=1:N

endfor;

figure

Answers.

0.5 1.0 1.5 2.0 2.5

0.2 0.4 0.6 0.8 1.0

figure

plot(ts,ys);

title ("y(t),y’(t), y(0)=0,y’(0)=1");

Answers.

0.5 1.0 1.5 2.0 2.5

0.2 0.4 0.6 0.8 1.0

Correct answer. A

(25)

Critical damping of a harmonic oscillator

Let

a= 2.5, b= 1, t0= 2.5, t1 = 2.5, (11)

ic₀= 0.8, ic₁ = 1.1, (12)

y⁰⁰(t) +ay⁰(t) +by(t) =f(t). (13)

y v

=A y

v

+B(f(t)) =

0 1

−b −a y v

+

0 1

(f(t)), (14)

~ y(0) =

y(0) v(0)

= ic₀

ic1

. (15)

Mathematica. {a,b, t0,t1, ic0,ic1}={?,?, ?,?, ?,?}

de={y’’[t]+a y’[t]+b y[t]==f[t]}

ic={y[t0]==ic0,y’[t0]==ic1}

A={{0,1},{-b,-a}}

Octave. clear all;

%graphics_toolkit ("gnuplot");

pkg load symbolic;

global a=? b=? A=[0,1;-b, -a];

t0=?; t1=?; ic0=?; ic1=-?;

syms y(t);

ydot=A*yvec;

endfunction N=100;

36, A

Exercise. Let f(t) = 0, (2) be true. How much isy(t1) ?

~

U_t₁_,t₀( ic0

ic₁

).

(26)

Exercise. Let f(t) = 0, (2) be true. Let U_t₁_/6,t₀(~z) = ~y(t1/6), where ~y(t) is the solution of (4) with the initial condition~y(t0) =~z. ComputeU_t₁_/6,t₀(~z), if~z=n1(0.2,0)^T+n2(0,0.2)^T, n1,2∈ {0,1,2,3,4,5}. Plot these points!

Mathematica. m=N[MatrixExp[(t1/6)*A]]

Octave. figure

u=expm((t1/6)*A);

endfor endfor;

hold on;

hold off;

title("[y(t1),y’(t1)] -> [y(t2),y’(t2)]")

(27)

Answers.

How much isspring(t₁) ?

y=x(1); v=x(2); s=x(3); f=x(4);

res(1)=ydot-v;

res(2)=vdot+s+f;

res(3)=s-b*y;

res(4)=f-a*v;

endfunction N=100;

(28)

y(t) =C₁e^λt+C₂te^λt. How much isC1, if (2) is satisfied?

Mathematica. lam=Sqrt[b];

ans2=ic1 - lam ic0;

Octave. ic1-sqrt(b)*ic0

Answers. A: 0.3 Correct answer. A 41, Ba

~

y(t) =C₁e^λ¹^t~v₁+C₂e^λ²^t~v₂, How much isC1(~v2)1, if (10) is satisfied?

Mathematica. jd=JordanDecomposition[A]

S=jd[[1]]

J=jd[[2]]

ans2a=Coefficient[

((S.MatrixExp[t J].Inverse[S]).

{{ic0},{ic1}})[[1]],t][[1]]/.t->0

(* OR (this might yield different result, due to the instability of the Jordan decomposition)*)

Coefficient[((MatrixExp[t A]).{{ic0},{ic1}})[[1]],t][[1]]/.t->0 (* OR *)

lam=Sqrt[b]

ans2=ic1-lam ic0

Octave. ic1-sqrt(b)*ic0

Answers. A: 0.3 Correct answer. A 42, D

(29)

figure

plot(ts,ys(:,1));

title ("y(t)");

Answers. ^0.5 ^1.0 ^1.5 ^2.0 ^2.5

0.5 0.6 0.7 0.8 0.9 1.0

figure

plot(ts,ys(:,1));

Answers.

0.5 1.0 1.5 2.0 2.5

0.5 1.0

Correct answer. A

(30)

44, F

figure

plot(ys(:,1),ys(:,2));

title ("t->[y(t),y’(t)]");

Answers.

0.6 0.7 0.8 0.9 1.0 0.5

1.0

γ : [t₀, t₁]→R³, γ(t) = (t, y(t), v(t))^T parametric curve!

(31)

figure

plot3(ts,ys(:,1),ys(:,2));

title ("t->[t,y(t),y’(t)]");

Answers.

0

1

2

0.6 0.8

1.0 0.0

0.5 1.0

figure

title ("[y;v] -> A[y;v]");

(32)

Answers. ^-^1.0 ^-^0.5 ^0.0 ^0.5 ^1.0

-1.0 -0.5 0.0 0.5 1.0

Show[vPlot,cPlot]

plot(ys(:,1),ys(:,2));

hold off;

(33)

Answers. ^-^1.0 ^-^0.5 ^0.0 ^0.5 ^1.0

-1.0 -0.5 0.0 0.5 1.0

ans9[[1,1]]

Plot[{U[[1,1]],U[[2,1]]},{t,t0,t1}]

for i=1:N

endfor;

figure

(34)

Answers.

0.5 1.0 1.5 2.0 2.5

-0.4 -0.2 0.2 0.4 0.6 0.8 1.0

figure

plot(ts,ys);

title ("y(t),y’(t), y(0)=1,y’(0)=0");

Answers.

0.5 1.0 1.5 2.0 2.5

-0.4 -0.2 0.2 0.4 0.6 0.8 1.0

Plot[{U[[1,2]],U[[2,2]]},{t,t0,t1}]

(35)

for i=1:N

endfor;

figure

Answers.

0.5 1.0 1.5 2.0 2.5

0.2 0.4 0.6 0.8 1.0

figure

plot(ts,ys);

title ("y(t),y’(t), y(0)=0,y’(0)=1");

Answers.

0.5 1.0 1.5 2.0 2.5

0.2 0.4 0.6 0.8 1.0

Correct answer. A

(36)

Distributions (δ, δ

⁰

)

Letχ_A(t) = 1, ift∈A, otherwise χ_A(t) is zero.

n= 1, a= 5, φ(t) =e^at, αn(t) =nχ[−1/(2n),1/(2n)](t), βn(t) = 1

√2πσn

e^−t²^/(2σⁿ²⁾, σn= 1/(4n),

˜

α_n(t) =n²(χ_[−1/n,0](t)−χ_[0,1/n](t)), β˜_n(t) =β_n⁰(t), hf, φi=

Z ∞

−∞

f(t)φ(t)dt

Mathematica. alpha[n_,t_]:=n HeavisideTheta[t+1/(2n)]HeavisideTheta[-t+1/(2n)]

alphat[n_,t_]:=n^2(-HeavisideTheta[-t+1/n]HeavisideTheta[t]+

HeavisideTheta[t+1/n]HeavisideTheta[-t]) sigma[n_]:=1/(4n)

beta[n_,t_]:=1/(sigma[n]Sqrt[2Pi])Exp[-0.5(t/sigma[n])^2]

betat[n_,t_]:=D[beta[n,s],s]/.s->t Octave.

53, A

Exercise. Plot αn(t) and βn(t) !

Mathematica. Plot[{alpha[n,t],beta[n,t]},{t,-1,1},PlotRange->All]

Octave.

Answers. ^-^1.0 ^-^0.5 ^0.5 ^1.0

2 4 6 8

Exercise. Plot βk(t), k=n−1, n, n+ 4 !

Mathematica. Plot[{beta[n-1,t],beta[n,t],beta[n+4,t]},{t,-1,1},PlotRange->All]

Octave.

(37)

Answers. ^-^0.3 ^-^0.2 ^-^0.1 ^0.1 ^0.2 ^0.3

2 4 6 8 10 12 14

Correct answer. A 55, B

Exercise. Plot αe_n(t) and βe_n(t)!

Mathematica. Plot[{alphat[n,t],betat[n,t]},{t,-1,1},PlotRange->All]

Octave.

Answers.

-1.0 -0.5 0.5 1.0

-100 -50 50 100

Correct answer. A 56, C

Exercise. Compute φ(0)− hα, φi ?

Mathematica. Exp[a*0]-NIntegrate[Exp[a t]n,{t,-1/(2 n),1/(2n)}]

Octave.

Answers. A: -0.0016675 Correct answer. A

(38)

57, D

Exercise. Compute (−φ⁰(0))− hα,φie ?

Mathematica. -(a) -(-NIntegrate[Exp[a t]n^2,{t,0,1/n}]

+NIntegrate[Exp[a t]n^2,{t,-1/ n,0}])

Octave.

Answers. A: 0.00333778 Correct answer. A 58, E

Exercise. Compute φ(0)− hβ, φi !

Mathematica. 1-NIntegrate[Exp[a t]beta[n,t],{t,-10,10}]

Octave.

Answers. A: -0.00125078 Correct answer. A 59, F

Exercise. Compute (−φ(0))− hβ, φie !

Mathematica. -a-NIntegrate[Exp[a t]betat[n,t],{t,-10,10}]

Octave.

Answers. A: 0.00125078 Correct answer. A 60, G

Exercise. Let

n= 2ⁱ, i= 0, . . . ,4, errn=|φ(0)− hα_n, φi|.

Plot the points [ln(n),ln(errn)] and find the best fitting line by linear regression (i.e. by the method of least squares)!

Mathematica. err=Table[{N[Log[2^i]],

Log[Abs[1-NIntegrate[Exp[a t]2î,{t,-1/(2 2î),1/(2 2î)}]]]},{i,0,4}]

model=LinearModelFit[err,t,t]

Show[ListPlot[err], Plot[model["BestFit"], {t, 0, 10}]]

Octave.

(39)

Answers.

0.5 1.0 1.5 2.0 2.5

-8 -6 -4 -2

Exercise. Let

n= 2ⁱ, i= 0, . . . ,4, errn=|φ(0)− hα_n, φi|.

Take the points [ln(n),ln(errn)] and find the best fitting line by linear regression (i.e. by the method of least squares)! What is the slope of that line?

Log[Abs[1-NIntegrate[Exp[a t]2î,{t,-1/(2*2î),1/(2*2î)}]]]},{i,0,4}]

Octave.

Answers. A: -2.00401 Correct answer. A 62, Hb

Exercise. Let

n= 2ⁱ, i= 0, . . . ,4, errn=|φ(0)− hγ_n, φi|, where

γ_n(t) =χ_[0,1/n](t)

Take the points [ln(n),ln(err_n)] and find the best fitting line by linear regression (i.e. by the method of least squares)! What is the slope of that line?

Log[Abs[1-NIntegrate[Exp[a t]2^i,{t,0,1/( 2^i)}]]]},{i,0,4}]

Octave.

Answers. A: -1.11748 Correct answer. A

(40)

63, I

Exercise. Let

n= 2ⁱ, i= 0, . . . ,4, rrn=|(−φ⁰(0))− hαen, φi|.

Plot the points [ln(n),ln(errn)] and find the best fitting line by linear regression (i.e. by the method of least squares)!

Mathematica.

Octave.

Answers.

0.5 1.0 1.5 2.0 2.5

-8 -6 -4 -2

Exercise. Let

n= 2ⁱ, i= 2, . . . ,13, err_n=|φ(0)− hα_n, φi|.

Take the points [ln(n),ln(err_n)] and find the best fitting line by linear regression (i.e. by the method of least squares)! What is the slope of that line?

Mathematica. err=Table[{N[Log[2^i]],Log[Abs[

-a-(-NIntegrate[Exp[a t](2^i)^2,{t,0,1/2^i}]

+NIntegrate[Exp[a t](2^i)^2,{t,-1/ 2^i,0}])]]},{i,0,4}]

ans=Coefficient[model["BestFit"],t]

Octave.

Answers. A: -2.01072 Correct answer. A 65, K

Exercise. y⁰(t) =δ(t), y(−1) = 0, y⁰_n(t) =αn(t), yn(−1) = 0.Plot y, yn-t!

Mathematica. ysol=DSolve[{y[-1]==0,y’[t]==DiracDelta[t]},y,t]

yf=ysol[[1,1,2,2]]

ynsol=DSolve[{y[-1]==0,y’[t]==alpha[n,t]},y,t]

ynf=ysol[[1,1,2,2]]

Plot[{yf,ynf},{t,-0.5,0.5}]

(41)

Octave.

Answers. ^-^0.4 ^-^0.2 ^0.2 ^0.4

0.2 0.4 0.6 0.8 1.0

Exercise. y⁰(t) =δ(t), y(−1) = 0, y⁰_n(t) =β_n(t), y_n(−1) = 0.Plot y, y_n-t!

Mathematica. ysol=DSolve[{y[-1]==0,y’[t]==DiracDelta[t]},y,t]

yf=ysol[[1,1,2,2]]

ynsol=DSolve[{y[-1]==0,y’[t]==beta[n,t]},y,t]

ynf=ysol[[1,1,2,2]]

Plot[{yf,ynf},{t,-0.5,0.5}]

Octave.

Answers. ^-^0.4 ^-^0.2 ^0.2 ^0.4

0.2 0.4 0.6 0.8 1.0

Exercise. y⁰⁰(t) =δ(t), y(−1) = 0, y⁰(−1) = 0, y⁰⁰_n(t) =α_n(t), y_n(−1) = 0, y⁰(−1) = 0.Plot y, y_n-t!

Mathematica. ysol=DSolve[{y[-1]==0,y’[-1]==0,y’’[t]==DiracDelta[t]},y,t]

yf=ysol[[1,1,2,2]]

(42)

ynsol=DSolve[{y[-1]==0,y’[-1]==0,y’’[t]==alpha[n,t]},y,t]

ynf=ysol[[1,1,2,2]]

Plot[{yf,ynf},{t,-0.5,0.5}]

Octave.

Answers. ^-^0.2 ^-^0.1 ^0.1 ^0.2

0.05 0.10 0.15 0.20

Correct answer. A 68, N

Exercise. y⁰⁰(t) =δ(t), y(−1) = 0, y⁰(−1) = 0, y⁰⁰_n(t) =β_n(t), y_n(−1) = 0, y⁰(−1) = 0.Plot y, y_n-t!

yf=ysol[[1,1,2,2]]

ynsol=DSolve[{y[-1]==0,y’[-1]==0,y’’[t]==beta[n,t]},y,t]

ynf=ysol[[1,1,2,2]]

Plot[{yf,ynf},{t,-0.5,0.5}]

Octave.

Answers. ^-^0.2 ^-^0.1 ^0.1 ^0.2

0.05 0.10 0.15 0.20

Correct answer. A

(43)

69, O

Exercise. y⁰⁰(t) =δ⁰(t), y(−1) = 0, y⁰(−1) = 0, y⁰⁰_n(t) =βe_n(t), y_n(−1) = 0, y⁰(−1) = 0.Plot y, y_n-t!

yf=ysol[[1,1,2,2]]

ynsol=DSolve[{y[-1]==0,y’[-1]==0,y’’[t]==beta[n,t]},y,t]

ynf=ysol[[1,1,2,2]]

Plot[{yf,ynf},{t,-0.5,0.5}]

Octave.

Answers. ^-^0.2 ^-^0.1 ^0.1 ^0.2

0.2 0.4 0.6 0.8 1.0

Correct answer. A 70, P

Exercise. y⁰⁰(t) = 0.4δ(t) + 1.6δ⁰(t), y(−1) = 0, y⁰(−1) = 0, y_n⁰⁰(t) = 0.4β(t) + 1.6βen(t), yn(−1) = 0, y⁰(−1) = 0.

Ploty, y_n-t!

Mathematica. ysol=DSolve[{y[-1]==0,y’[-1]==0,

y’’[t]==0.4DiracDelta[t]+1.6D[DiracDelta[t],t]},y,t]

yf=ysol[[1,1,2,2]]

ynsol=DSolve[{y[-1]==0,y’[-1]==0,

y’’[t]==0.4beta[n,t]+1.6betat[n,t]},y,t]

ynf=ysol[[1,1,2,2]]

Plot[{yf,ynf},{t,-0.5,0.5}]

Octave.

(44)

Answers. ^-^0.3 ^-^0.2 ^-^0.1 ^0.1 ^0.2 ^0.3

0.5 1.0 1.5

Correct answer. A 71, Q

Exercise. y⁰⁰(t) = −3y(t)−y⁰(t) +δ(t), y(−1) = 0, y⁰(−1) = 0, y_n⁰⁰(t) = −3y_n(t)−y_n⁰(t) +β_n(t), y_n(−1) = 0, y⁰_n(−1) = 0.Plot y, yn-t!

Mathematica. ysol=DSolve[{y[-1]==0,y’[-1]==0,

y’’[t]==0.4DiracDelta[t]+1.6D[DiracDelta[t],t]},y,t]

yf=ysol[[1,1,2,2]]

ynsol=DSolve[{y[-1]==0,y’[-1]==0,

ynf=ysol[[1,1,2,2]]

Plot[{yf,ynf},{t,-0.15,0.15}]

Octave.

Answers. ^-^0.15 ^-^0.10 ^-^0.05 ^0.05 ^0.10 ^0.15

0.02 0.04 0.06 0.08 0.10 0.12 0.14

Correct answer. A 72, R

Exercise. y⁰(t) =−3y(t) +δ(t), y(−1) = 0, y_n⁰(t) =−3y_n(t) +βn(t), yn(−1) = 0.Plot y, yn-t!

(45)

Mathematica. ysol=DSolve[{y[-1]==0, y’’[t]==DiracDelta[t]},y,t]

yf=ysol[[1,1,2,2]]

ynsol=DSolve[{y[-1]==0,y’[-1]==0,

ynf=ysol[[1,1,2,2]]

Plot[{yf,ynf},{t,-0.15,0.15}]

Octave.

Answers. ^-^0.2 ^-^0.1 ^0.1 ^0.2

0.2 0.4 0.6 0.8 1.0

Correct answer. A 73, S

Exercise. y⁰(t) =−3y(t) +δ(t), y(−1) = 0, y_n⁰(t) =−3y_n(t) +β_n(t), y_n(−1) = 0.Plot y, y_n-t!

Mathematica. ysol=DSolve[{y[-1]==0, y’’[t]==DiracDelta[t]},y,t]

yf=ysol[[1,1,2,2]]

ynsol=DSolve[{y[-1]==0,y’[-1]==0,

ynf=ysol[[1,1,2,2]]

Plot[{yf,ynf},{t,-0.15,0.15}]

Octave.

(46)

Answers. ^-^0.2 ^-^0.1 ^0.1 ^0.2

0.2 0.4 0.6 0.8 1.0

Correct answer. A

(47)

StudentNumber: 91

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

A A A A A A A A A A A A A A A A A A A A

21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40

41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80

A A A A A A A A A A A A A